Mapping Soil Erosion in a Quaternary Catchment in Eastern Cape Using Geographic Information System and Remote Sensing

In South Africa, soil erosion is considered as an environmental and social problem with serious financial implications particularly in some rural areas where this geomorphological phenomenon is widespread. An example is the Umzimvubu Local Municipality, where most households are strongly reliant on agriculture for their livelihood. Sustainable agriculture and proper land management in these rural areas require information relevant to the spatial distribution of soil erosion. This study was therefore aimed at generating such information using Landsat8 Operational Land Imager (OLI)derived vegetation indices (VIs) including the Normalised Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), as well as Soil and Atmospherically Resistance Vegetation Index (SARVI). Raster calculator in ArcMap10.2 was used to classify soil erosion features based on selected suitable thresholds in each VI. SPOT6/7 (Satellites Pour l’Obsevation de la Terre) multispectral data and Google Earth images were used for ground truth purposes. SAVI achieved the highest overall classification accuracy of 83% and kappa statistics of 64%, followed by NDVI and SARVI with equal overall accuracy of 81% and slightly different kappa statistics of 60% for the former and 59% for the latter. Using these indices, the study successfully mapped the spatial distribution of soil erosion within the study area albeit there were some challenges due to coarser spatial resolution (15mx15m) of Landsat8 image. Due to this setback, image fusion and pan-sharpening of Landsat8 with higher spatial resolution images is strongly suggested as an alternative to improve the Landsat8 spatial


Introduction
Soil erosion is considered one of the world's most critical environmental concerns owing to its adverse effects on both the natural environment and human society (Jie et al., 2002;Le Roux et al., 2008;Aiello et al., 2015).Globally, approximately 2 billion hectares (22 percent) of all cropland, pasture, forest and woodland around the world have been degraded since last century (Oldeman et al., 1991;Scherr and Yadav, 1997;Pimentel, 2006).In financial terms, according to Scherr (1999) cited in Jie et al. (2002), at least $28 billion have been lost per year due to soil degradation in drylands around the world.While soil erosion affects both the developed and developing world, the latter appears to have experienced more severe erosion than its developed counterpart.In South Africa, more than 70% of the land is subject to significant levels of soil erosion (Garland et al., 2000).Many soil erosion-borne studies (including Beckedhal and De Villiers, 2000;Boardman et al., 2003;Ngetar, 2011) have been conducted over the past few years in South Africa.However, one weakness of the South African soil erosion research is the limited information on where the problem is highly concentrated (Le Roux et al., 2007).Considerable attention has been paid to the "what" question rather than to the "where" question and this calls for more research to generate adequate information regarding the spatial distribution of the problem.
Studies making use of geospatial technologies including GIS and remote sensing to map the spatial distribution of soil erosion in South Africa are relatively recent (Le Roux and Sumner, 2012;Seutloali et al., 2015).For example, Mararakanye and Nethengwa (2012) mapped the spatial distribution of gully erosion in Limpopo using remote sensing and traditional techniques.Le Roux and Sumner (2012) also digitised continuous and discontinuous gullies from SPOT5 imagery showing the spatial extent of gullies at a catchment scale.More recently, using Landsat multispectral data, Seutloali et al. (2016) produced soil erosion severity index maps of the former homelands of Transkei.Despite the use of GIS and remote sensing in soil erosion studies in South Africa, there is still a need to map the spatial extent of the problem (Le Roux et al., 2007) in other erosion-prone areas such as the Umzimvubu Municipality (Eastern Cape) within which the study area (catchment) falls.
Dissemination of information relevant to the spatial distribution of soil erosion within the municipality is critical to land conservation and management.If governance is about decisionmaking, then up-to-date, accurate, complete and usable information on soil erosion is indispensable to governance (Macharia, 2005) at a local level, namely, the municipality.However, this requires the availability of high spatial resolution imagery.In the absence of such imagery, freely available coarse resolution imagery can be used to provide an insight into the nature of soil erosion in an area.
South African Journal of Geomatics, Vol. 6. No. 1, April 2017 13 Different techniques exist for the assessment of soil erosion phenomenon including the most widely used Universal Soil Loss Equation (USLE) developed by Wischmeier and Smith (1965).The use of USLE and its enhanced versions have been instrumental in estimating soil loss rates with reasonable costs and better accuracy especially when combined with GIS and remotely sensed data (Owusu, 2012;Alexakis et al., 2013;Ganasri and Ramash, 2015).Nonetheless, as highlighted by Merritt et al. (2003) and Kinnell (2010), USLE (including its enhanced versions) was designed to assess rill and inter-rill erosion and hence its application is problematic in areas where ephemeral and classical gullies are dominant forms.Besides this apparent limitation, USLE or any other soil loss model is not considered in this paper since the paper is primarily concerned with the spatial distribution of the erosion phenomenon as opposed to modelling absolute values of soil loss.Various remote sensing techniques (automated and semi-automated) exist for mapping soil erosion, ranging from pixel-based to object-oriented image classification with the latter relying on fine spatial resolution images.While the growing availability of high spatial resolution images such as IKONOs, QuickBird, and WorldView has facilitated a shift from traditional pixel-based to object-based methods (Shruthi et al., 2011;Mayr et al., 2016), these images are not readily or freely available.
Additionally, such high spatial resolution images have limited spectral bands as noted by Taruringa (2008) and Shruthi et al. (2011).It is therefore not surprising that medium resolution images such as Landsat still appeal to some soil erosion researchers (Seutloali et al., 2016;Dube et al., 2017).
The use of semi-automated image classification techniques such as vegetation indices (VIs) in particular have gained momentum in soil erosion research over the past few decades (Mathieu et al., 1997;Singh, 2004).Since its introduction in the 1970s by Rouse et al. (1973), the Normalised Difference Vegetation Index (NDVI) has been used quite extensively in soil erosion-related research (Vaidyanathan et al., 2002;Li et al., 2010;Seutloali et al., 2016).However, the extensive use of NDVI has also presented significant weaknesses.Govaerts and Verhulst (2010) note that satellitebased NDVI results are subject to interference by non-vegetation factors including, but not limited to atmospheric conditions and soil background.Consequently, various modifications have been proposed to address the sensitivity of NDVI to non-vegetation factors (Lawrence and Ripple, 1998).
The Soil Adjusted Vegetation Index (SAVI) proposed by Huete (1988), and Soil and Atmospherically Resistant Vegetation Index (SARVI) developed by Huete and Liu (1994) are amongst the most widely used modifications of NDVI in soil erosion research.
In this paper, satellite image-derived NDVI, SAVI, and SARVI are considered mainly because they are some of the simplest, cheapest, and quickest feature extraction techniques (Singh, 2004;South African Journal of Geomatics, Vol. 6. No. 1, April 2017 14 Gandhi et al., 2015;Alhawiti and Mitsova, 2016;Sonawane and Bhagat, 2017) that can be used to map soil erosion events.Soil erosion is a dynamic process requiring constant monitoring while keeping up-to-date information on its spatial distribution.However, in the study area and the case of Umzimvubu Municipality at large, existing information on soil erosion (Madikizela, 2000) is outdated and sparse.Apparently, there are no soil erosion studies in the area which have taken advantage of GIS and remote sensing and their assessment methods, particularly the use of VIs.Therefore, the main objective of this study is twofold: (1) to map the spatial distribution of soil erosion in the study area using three Landsat-derived VIs, namely; NDVI, SAVI, and SARVI, and (2) to assess the accuracy of VI-derived soil erosion maps and determine the best VI for detecting soil erosion features at a catchment level.

Description of the Study Area
The study area constitutes one of the most severely eroded catchments in the Eastern Cape and is located in the central part of Umzimvubu Local Municipality.It includes sections of the Matatiele and Ntabankulu local municipalities on the north and south, respectively.The approximate geographic coordinates of the catchment area are 30° 38´ 23.5´´-30° 55´ 29.50´´ S and 28° 42´ 36.65´´-29° 10´ 19.17´´ E (Figure 1), and covers approximately 503km² of the total municipal surface area (2 506km²).The study area exhibits a highly uneven topography with an elevation range of 820m to 1970m above sea level, for low-lying and high-lying areas, respectively.

Materials
The image data used in this study include Landsat8 Operational Land Imager (OLI) and Satellites Pour l'Obsevation de la Terre (SPOT6/7) multispectral images obtained in slightly different acquisition dates.The term "SPOT6/7" which will be used throughout this paper refers to a newly launched SPOT image product which combines the capabilities of SPOT6 and SPOT7 imagery.
Landsat8 OLI consists of nine multispectral bands with a spatial resolution of 30m for bands 1 to 7 and 9. Band 8 (panchromatic) has a spatial resolution of 15m whereas bands 10 and 11 (thermal bands) are collected at 100m spatial resolution.SPOT6/7 consists of four multispectral bands with a spatial resolution of 6m and a panchromatic band with a spatial resolution of 1.5m.Landasat8 OLI was acquired in 2014-07-20 at 07:57:13 whereas SPOT6/7 was acquired in 2014-06-09 at 07:49:31 according to their respective scene centre times.SPOT6/7 and Google Earth were used as ancillary Landsat8 data for ground truth purposes.Besides providing superior spatial resolution, Google Earth was chosen as a surrogate for field work, thereby avoiding travelling costs.SPOT6/7 and Landsat8 OLI were chosen because they are readily available free of charge, but most importantly because of their suitability for detecting erosion features (Taruvinga, 2008;Dube et al., 2017).Moreover, Seutloali et al. (2016) maintain that Landsat data has proven to be significant in examining soil erosion occurrence especially in instances where in-depth field work remains a challenging task.
Landsat8 OLI was freely downloaded from the United States Geological Survey (USGS) website, while SPOT6/7 data was freely obtained from the South African Space Agency (SANSA).The Red (band 4), and near-infrared (NIR) (band 5) of Landsat8 OLI were utilised to compute NDVI and SAVI.SARVI was generated from the Blue (band 2), Red, and NIR bands of Landsat8 OLI.In order to minimise the effects of soil variability and atmospheric noise on vegetation, a soil calibration factor of 0.5 was used for both SAVI and SARVI including a blue-band normalisation for the latter (Huete et al., 1988;Huete and Liu, 1994).

Pre-processing
Prior to the actual image processing, Landsat8 OLI and SPOT6/7 scenes were subset to the study area (catchment) in order to restrict analysis; thereby reducing time required for image processing (Campbell and Wynne, 2011).Both SPOT6/7 and Landsat8 OLI multispectral images were obtained already pre-processed by suppliers.However, in order to obtain maximum surface reflectance, Landsat8 digital numbers (DN) were converted to Top Of Atmosphere (TOA) reflectance using the radiometric rescaling coefficients provided in the Landsat8 product metadata file (MTL file).All calculations were performed using the "raster calculator" tool within the ArcMap 10.2 environment.
The following equation [1] available in the USGS Landsat8 website was used to convert digital numbers (DN) to TOA reflectance (Alhawiti and Mitsova, 2016): Where ΄ is the TOA planetary reflectance (without correction for sola angle),   is the Bandspecific multiplicative rescaling factor,   is the Band-specific additive rescaling factor, and   is the Quantized and calibrated standard product pixel values (DN).
Three Landsat8 OLI multispectral bands relevant to the computation of VIs, viz.Band 2 (Blue), Band 4 (Red), and Band 5 (NIR) were pan-sharpened to 15m spatial resolution using the panchromatic image (Band 8) of Landsat8 OLI.While various pan-sharpening algorithms have been proposed in the literature, the most widely used include Resolution Merge, Modified IHS (intensity hue saturation) Resolution Merge, and Adaptive IHS Pan-sharpening, amongst others (Rahmani et al., 2010;Zhang and Mishra, 2012).In this paper, the Resolution Merge pan-sharpening was considered based on the obtained output results relative to other algorithms that were tried and tested.
Pan-sharpening was conducted in ERDAS IMAGINE 2013 software using the built-in "Resolution Merge" module.

VIs Computation
The red and near-infrared (NIR) bands are probably the most important bands featured in all VIs calculations.This is so because vegetation reflects high in the near-infrared and low in the red band.
By maximising this linear relationship between the red and NIR bands, VIs allow for sharp spectral discrimination between vegetation and non-vegetation attributes including soil erosion (Taruvinga, 2008).Given this linear relationship where soil erosion increases with decreasing vegetation (Seutloali and Beckedahl, 2015), an assumption can be made that areas lacking vegetation cover represents erosion, and would therefore exhibit a particular VI value or range (threshold) that could be used to identify eroded areas (Taruvinga, 2008).Furthermore, once erosion features are identified through a threshold technique, such erosion features can be mapped (Cyr et al., 1995) Where   is the near-infrared radiation,   is the visible red radiation, and L is the soil adjustment factor.An L value of 0.5 minimises the effects of soil brightness variations and eliminates the need for additional calibration for different soils (Huete and Liu, 1994).
Where   is   − (  −   ),   is the near-infrared radiation,   is the visible red radiation, and   is the visible blue radiation, L is the soil adjustment factor.Huete and Liu (1994) combined the soil adjustment factor (L) and a blue-band normalisation to correct for both soil and atmospheric noise.

Soil Erosion Classification.
Three VIs including NDVI, SAVI, and SARVI were produced mainly for soil erosion classification within the study area.All VIs were generated in ERDAS IMAGINE 2013 using predetermined formulas.The "arbitrary profiler" (spatial profiler) tool in ENVI5.2 was used to extract soil erosion reflectance values from selected gullies within the study area.Using a procedure similar to that of Taruvinga (2008) and Gandhi et al. (2015), best thresholds for soil erosion classification were determined for each VI on a trial and error basis.The benchmark used to select the best threshold was determined based on the accuracy of the derived erosion maps (i.e.overall accuracy and kappa statistics results).Raster calculator in ArcMap10.2was then used to classify eroded areas based on the selected classification thresholds for each VI.

Accuracy Assessment
One of the simplest and yet widely accepted ways of assessing accuracy is to compare ground reference test data with remote sensing classification maps (Jensen, 2005).Integral to this, is obtaining the appropriate number of unbiased ground reference points.In this study, 100 random sample points were created using a "Create Random Points" tool in ArcMap 10.2.Within the "Create Random Points" tool dialogue box, the "Constraining Feature Class" option was set to the polygon shapefile of the study area in order to restrict random points within the boundaries of the study area.
It is very important to collect ground reference test information as close to the date of image acquisition as possible (Jensen, 2005;Campbell and Wynne, 2011).Accordingly, high resolution and most recent images of SPOT6/7 (2014) and Google Earth image (2013) were used for ground truth to assess the accuracy of the 2014 Landsat classified images.The GIS-derived random points were initially converted to kml (keyhole makeup language) and exported to Google Earth for ground truth purposes.Given that soil erosion features are very small and often linear in shape, most of the random points were widely distributed on non-erosion features, with very few sitting on erosion features.
Therefore, each random point was assigned to the adjacent erosion feature on Google Earth and these totalled to 40 while the remaining 60 random points sat on non-erosion features.
Confusion matrix was tabulated for each VI.Four levels of accuracy including user's accuracy, producer's accuracy, overall accuracy, and kappa statistics were calculated.Confusion matrix was chosen because it is one of the most widely used and accepted ways of expressing classification accuracy (Story and Congalton, 1986;Lillesand and Kiefer, 2000;Sonawane and Bhagat, 2017).The confusion matrix helped to examine the relationship between known reference data and the corresponding results of VI-derived maps.

Soil Erosion Classification from VIs
After a series of trial runs, 0.12 to 0.22 reflectance values were selected as the best threshold to classify soil erosion in NDVI while 0.15 to 0.25 and 0.13 to 0.22 were chosen for SAVI and SARVI, respectively.The "Raster Calculator" tool in ArcMap10.2was used to classify erosion features based on the above-mentioned thresholds.Shown in Figure 3 are classification results of Landsat8 images based on the selected threshold reflectance values for each VI.The distribution of soil erosion is overlaid on SPOT6/7 imagery (Figure 4).

Accuracy Assessment Results of VI-derived Maps
A necessary step after remote sensing image classification is the accuracy assessment of the results.
Table 1 provides the accuracy assessment report of VI-derived classification results with respect to producer's and user's accuracy along with overall accuracy and kappa statistics for erosion and nonerosion classes.Figure 5 provides a summary of the overall accuracy and kappa statistics.
The results indicate that NDVI achieved highest accuracy for producer's accuracy (78.4%) in the erosion class (Table 1).SAVI outperformed other VIs with an overall accuracy of 83% and kappa statistics of 0.64 (64%).NDVI and SARVI yielded the same results in terms of overall accuracy with slightly different kappa statistical results viz.0.60 (60%) for the former and 0.59 (59%) for the latter.
Relative to other VIs, SARVI recorded the lowest producer's accuracy of 62.5% and a commanding 86.7% user's accuracy for the erosion category.Since SAVI achieved the highest overall accuracy and kappa statistic results (Figure 5), it was chosen to classify soil erosion within the study area.
Figure 6 shows the distribution of soil erosion over the entire study area while Figure 7 is zoomed in to the western portion of the study area.

Discussion
One objective of this study was to map the spatial distribution of soil erosion using three Landsatderived VIs.Notwithstanding the success of the study in mapping the spatial distribution of soil erosion, the results ought to be examined from its proper perspective.Despite an attempt to improve the spatial resolution, not all soil erosion features have been classified, for example rills, sheet and some small gullies have been excluded possibly because their sizes fall below the spatial resolution (15m) of the pan-sharpened Landsat8 OLI.Within the study area, the spectral reflectivity of soil erosion varies considerably, and in some cases tends to be similar to non-erosion features (for example bare soil).The possibility therefore exists that some soil erosion features may have been classified as non-erosion features.For example, the spectral reflectance of sheet erosion resembles that of bare soil surfaces making it difficult if not impossible to spectrally discriminate between the two features.Admittedly, an overall accuracy of 83% for SAVI and 81% for NDVI and SARVI denotes a good or perhaps high level of accuracy (Taruvinga, 2008;Mararakanye and Nethengwe, 2012;Sonawane and Bhagat, 2017).However, these may be misleading as well, given that much of the study area is covered by the non-erosion features, for instance; at least 80% of the study area in each VI, constitutes non-erosion features.Accordingly, using an overall accuracy as a proxy of accuracy could be very biased because of the larger spatial coverage of the non-erosion class (Taruvinga, 2008).

Conclusion and Recommendations `
Overall, this study successfully mapped the spatial distribution of soil erosion in a quaternary catchment within the Umzimvubu Municipality in Eastern Cape, South Africa.The kappa statistical results of 0.64 for SAVI, 0.60 for NDVI, as well as 0.59 for SARVI which are all considered acceptable and valid, justify this conclusion.Relative to other indices, it has been proven that SAVI is the most suitable VI for mapping soil erosion at a catchment level.Despite these positive results, the detection of some small erosion features proved to be a challenge in this study partly because of the coarser spatial resolution (15mx15m) of Landsat8 OLI.This challenge was further compounded by the spectral homogeneity of soil erosion and non-erosion features.It is therefore highly recommended that soil erosion studies use alternative methods related to pan-sharpening and fusion of Landsat8 with panchromatic or multispectral data of higher spatial resolution images such as SPOT, QuickBird, IKONOS, and WorldView images amongst others.One way to improve the detection and spectral discrimination of erosion features is to integrate spectral and textural information through object-based image classification.Previous studies have shown that land managers and policy makers are more interested in the spatial distribution of soil erosion than in absolute values of soil loss (Lu et al., 2004).Therefore, this study sheds light not only on the spatial distribution of soil erosion in the study area but also provides a useful framework for land managers and decision makers when planning soil erosion management at the catchment level.

Acknowledgement
The authors would like to thank the United States Geological Survey (USGS) and the South African Space Agency (SANSA) for providing free Landsat and SPOT data respectively.

Figure 1 .
Figure 1.Location of the study area within the Umzimvubu Local Municipality.

Figure 7 .
Figure 7. Soil erosion distribution in the western section of the study area.
Where   is the near-infrared radiation, and   is the visible red radiation.
Symeonakis and Drake (2004)ies that have applied threshold technique include that ofMathieu et al. (1997),Vaidyanathan et al. (2002), andSymeonakis and Drake (2004), amongst others.In this paper, a threshold technique is utilised in the subsequent section to classify soil erosion with the aid of each VI.While all VIs namely; NDVI, SAVI, and SARVI generally rely on red and NIR bands to calculate vegetation, they are computed using different formulas as indicated in equations [2],[3], and [4].

Table 1 .
Accuracy assessment report for erosion and non-erosion classes in VIs.Figure 5. Summary statistical results of overall accuracy and kappa statistics.

Table 1
results, particularly kappa statistics, it is apparent that SAVI is the best VI for detecting soil erosion features at a catchment scale and a spatial resolution of 15m.While the kappa statistics of 0.64 (64%)