Landslide manual and automated inventories, and susceptibility mapping using LIDAR in the forested mountains of Guerrero, Mexico

ABSTRACT Landslides are a pervasive natural disaster, resulting in severe social, environmental and economic impacts worldwide. The tropical, mountainous landscape in South-West Mexico is predisposed to landslides because of frequent hurricanes and earthquakes. The main goal of this study is to compare landslide susceptibility maps in Guerrero derived using high-resolution LIDAR (light detection and ranging) data from both a manual landslide event inventory and an automated landslide inventorying algorithm. The paper also highlights the importance of applying LIDAR data in landslide inventorying and susceptibility mapping. We mapped landslides based on two approaches: (1) manual mapping using satellite images and (2) automatic identification of landslide morphology employing the Contour Connection Method (CCM). We produced a landslide susceptibility map by computing the probability of landslide occurrence from statistical relationships of inventoried landslides detected with LIDAR digital terrain models (DTMs) and derived landslide-causing factors using the logistic regression method. Our results suggest that the automated inventory derived through the CCM algorithm with LIDAR DTMs effectively minimizes the time-consuming and subjective manual inventorying process. The high overall prediction accuracy (up to 0.83) from logistic regression demonstrates the validity and applicability deriving reliable landslide susceptibility maps from an automated inventory; however, LIDAR data are required.


Introduction
Landslides are a pervasive natural disaster, imparting severe social, economic and environmental impacts worldwide. These hazards are often paired with extreme weather events and/or earthquakes, producing mass movements of rock and soil on otherwise metastable slopes, sometimes occurring rapidly with little warning. Landslides can be devastating for communities, costing lives and hurting economies, often through the destruction of buildings and infrastructure, disruption of communications, damage of gas lines, water and sewage, and destruction of natural and cultural heritage monuments. These phenomena are often triggered by extrinsic factors (such as extraordinary rainfall or earthquakes) and intrinsic factors (such as geologic conditions, slope, vegetation and human activity). Many landslides occur along river systems and provide significant volumes of material to the lowlands (foothills and plains). These colluvial and fluvial materials increase the destructive power of debris and mudflows along the river systems, increasing the risk for human settlements and activities (Atkinson & Massari 1998;Dai & Lee 2002;Ohlmacher & Davis 2003).
The slopes in the mountainous areas of Sierra Madre del Sur in Guerrero, South-West Mexico, are frequently subject to intense rain from hurricanes (e.g. Hurricane Manuel on September 2013 triggered numerous landslides) as well as strong shaking from earthquakes, resulting in high landslide activity. Thus, identification of areas susceptible to landslides is fundamental in land management and planning for human occupation in this region. However, despite previous research (e.g. Legorreta-Paulin et al. 2012;Ramirez-Herrera et al. 2012;Legorreta-Paulin et al. 2013, and many others) and government efforts (e.g. Mendoza L opez & Dominguez Morales 2004; Programa Nacional de Protecci on Civil 2014-2018) to quantify and monitor landslides and other processes related to landslides, areas highly susceptible to landslides have been difficult to assess accurately. The low spatial and temporal resolution of remote sensing data cannot adequately capture the high spatial and temporal variability of landslide occurrence at the level of details required. The enhanced accuracy and detail of LIDAR (Light Detection and Ranging) data are necessary to evaluate the development of potential landslides, debris and mudflows to establish a practical methodology of landslide inventorying and susceptibility mapping, together with hazard forecasting. Worldwide examples have shown the utility of topographic data derived from airborne or terrestrial laser scanning in landslide investigations (e.g. Chigira et al. 2004;Schulz 2007;Booth et al. 2009;Guzzetti et al. 2012;H€ olbling et al. 2012;Jaboyedoff, Choffet, et al. 2012;Jaboyedoff, Oppikofer, 2012;Pradhan et al. 2012;Van Den Eeckhaut & Herv as 2012;Chen et al. 2013;Wang et al. 2013;Chen et al. 2014;Jebur et al. 2014;Lin et al. 2014;Scaioni et al. 2014;Chen et al. 2015;Dou et al. 2015;Li et al. 2015;Mahalingam & Olsen 2015;Mahalingam et al. 2016). Note that these studies tend to focus either on landslide detection or on landslide susceptibility mapping based on an already produced landslide inventory, without combining or examining the linkages between those two stages.
Hence, the primary goal of this study is to explore the differences of automated landslide inventorying techniques with high-resolution LIDAR digital terrain models (DTMs) and manual inventories using post-hurricane satellite data in Guerrero to produce high-quality automated susceptibility maps. To examine the feasibility of this process, we selected two study regions in the Coyuca River basin within the mountains of Guerrero, covering an area of >20 km 2 ( Figure 1). First, a manual inventory was performed on landslides triggered by Hurricane Manuel using satellite imagery obtained after the event. To avoid subjectivity, we also identified landslide features automatically employing the Contour Connection Method (CCM; Leshchinsky et al. 2015). A landslide susceptibility map was then developed by computing the probability of landslide occurrence from statistical relationships of existing landslides with LIDAR-elevation-derived landslide-causing factors using logistic regression (Regmi et al. 2014). The applied approach enabled the identification of ca. 90% of the landslides within the high and very high susceptibility zone, suggesting that the method was applicable for evaluating landslide susceptibility in tropical mountains of South-West Mexico and other similar areas. Finally, we demonstrate how the inventory source affects the output susceptibility maps by comparing susceptibility models based on the two different inventories. Ultimately, we show that an automated landslide inventory can be immediately integrated into susceptibility mapping, automating the entire process; however, high-resolution LIDAR topographic data are required to obtain models with the highest precision and accuracy.

Study area
Throughout the Coyuca basin (Figure 1), loosely deposited, quaternary alluvial sediments lie along the main rivers. These deposits are prone to frequent landslides, intense rainstorms produced by hurricanes, seismic activity and anthropogenic land modification. The combination of these factors highlights the importance of evaluating this region for landslide susceptibility.
The Coyuca basin is located within the Sierra Madre del Sur mountain range, which runs parallel to the Pacific coast in South-West Mexico (Figure 1). The basin is situated along an active subduction zone where the Cocos and North America plates converge, generating earthquakes in coastal regions of the Guerrero state. Nonetheless, this region has experienced no significant thrust earthquakes since 1911 (Anderson et al. 1989;Kostoglodov & Ponce 1994). Despite the relative calm in recent history, it is plausible that a major inter-plate earthquake with magnitude Mw 8.1-8.4 may occur along this segment of the coast in coming years (Su arez et al. 1990). Precambrian to Cretaceous gneiss of the Xolapa Complex, together with Eocene to Miocene granite and granodiorite, dominate the area (Carranza et al. 1999).
The two selected study regions are situated in the northern part of the Coyuca River basin, covering an area of approximately 22.3 km 2 . The larger northern region (No. 1; 15.6 km 2 ) encompasses the upstream section of the Coyuca River (Figure 1), characterized by rugged topography with relatively high relief and a dendritic drainage pattern. The elevation ranges from ca. 820 m, to more than 1530 m above mean sea level (amsl). Within the northern part of this region, the village of La Pintada witnessed a large deep-seated landslide triggered by Hurricane Manuel on 16 September 2013, resulting in 71 fatalities. The smaller region No. 2 (6.7 km 2 ) includes the E-W Coyuca River valley, with steep south-facing slopes and gentler north-facing slopes. The elevation varies here from »320 m (valley floor) to about 650 m amsl at the peaks.
Slopes in both study regions are covered with tropical and subtropical forests, most of which grow in poorly mixed, shallow gravelly or clayey soil deposits with little or no profile development. The northern part of region No. 1 near La Pintada village and the slopes facing south in region No. 2 are occupied mainly by pastureland and agriculture (Figure 1).
The climate within the region is classified as humid to sub-humid, warm to very warm (IG-UNAM 2007), with annual precipitation average rates ranging from ca. 600 mm/yr close to the town of Coyuca, and up to >1100 mm/yr in the central mountainous part of the Coyuca basin, closer to the shore ( Figure 2). Generally, higher rates of precipitation are observed in higher altitudes. Rainfall typically occurs during a five-month rainy period in the region, often during June to October, with the maximum rates in September reaching ca. 300 mm (Figure 2). The higher precipitation averages in September are often associated with the occurrence of hurricanes. Hurricane Manuel, a Category 1 hurricane (on the Saffir-Simpson Hurricane Wind Scale) made landfall as a tropical storm on the coast of Guerrero from the 13 to 19 September 2013, and resulted in 123 deaths in Mexico (Pasch & Zelinsky 2014) and numerous landslides. The total amount of rain recorded in Acapulco during this event reached more than 1100 mm (Pasch & Zelinsky 2014), which is approximately equal to the average annual precipitation for the entire Coyuca basin ( Figure 2).

Methods and materials
In this study, we produced landslide inventories using two different approaches including (1) a manual identification of landslides using satellite data collected after Hurricane Manuel (event landslide inventory), and (2) an automated identification of landslide morphology employing CCM (historical landslide inventory) with LIDAR-derived DTMs. We analysed landslide susceptibility using logistic regression utilizing each inventory separately. Data-sets to describe the landslide-causing factors were derived from LIDAR DTMs.
The LIDAR topographic data was collected on 19 March 2015 using a CESSNA TU206H aircraft with a Riegl Q-780 flying at the elevation of 700 m above ground at a nominal density of eight points per metres and altimetric accuracy of 35 cm. The Riegl Q-790 has a laser pulse repetition rate of 400 Hz and field of view of 60 (+30 /¡30 ). Riegl Ò RiProcess (http://www.riegl.com/products/soft ware-packages/) was used to produce a geo-referenced point cloud. The resulting point cloud was classified into ground and other basic categories (e.g. vegetation, infrastructure, buildings, etc.) using TerraScan Ò and TerraModeler Ò (https://www.terrasolid.com/home.php), followed by manual verification. A Triangulated Irregular Network (TIN) was generated for the bare earth ground classified point cloud data. For improved computational efficiency, a 1-m resolution raster version of the DTM was created from sampling the TIN.

Landslide inventory
Landslide inventorying is a critical initial step to enhance susceptibility mapping because it highlights the distribution, types and patterns of past landslides (Brabb 1991;Malamud et al. 2004;Burns & Madin 2009;Guzzetti et al. 2012). The delineation of the spatial extents of past landslides has been improved significantly in past years by the increasing availability of LIDAR, primarily because many of the morphological features associated with landslides are more apparent (e.g. Lin et al. 2014;Chen et al. 2015;Grejner-Brzezinska et al. 2015). Numerous methods to identify landslides using LIDAR topographic data have been proposed (e.g. Guzzetti et al. 2012;Scaioni et al. 2014). Despite these improvements, the landslide inventorying process is still challenging as manual digitizing of landslide extents can be time consuming and subjective (Burns & Madin 2009). Computer automation can expedite this process with more objectivity, particularly by replicating the manual digitizing process, namely by searching for headscarp morphology and hummocky topography associated with landslide deposits (Schulz 2004;Burns & Madin 2009). In this study, we produced landslide inventories employing two approaches, including manual mapping using satellite images, and automatic identification of landslide morphology employing CCM algorithm (Leshchinsky et al. 2015).
Characterization of general statistics describing landslide shape was performed for identified landslides following a modified procedure proposed by Niculiţǎ (2016). A Bounding Box tool in Arc-Map 10.2.2 Ò (http://desktop.arcgis.com/en/) was adapted to use azimuthal congruence with the LIDAR-derived aspect (slope direction) raster. Width (W) and length dimensions of each landslide (L) were defined as the dimensions perpendicular and parallel to the downslope direction of the landslide, respectively. Although this may demonstrate somewhat skewed values when landslides deposits extend through channels that are tortuous, it provides a reasonable means of evaluating generalized landslide dimensions as an aggregate.
For the accuracy assessment of both inventories, manual and automatic, we selected a 1 km 2 test section denoted by a polygon. In this area, we manually identified landslide extends using the LIDAR DTM, which were subsequently verified in the field. This manual image was treated as ground truth data (i.e. reference data). For the validation of manual inventory, we used only landslides produced by the Hurricane Manuel, whereas for the automatic inventory using CCM, we utilized all identified features. The accuracy assessment consisted of a comparison between reference data and both inventories. Apart from a visual comparison, we calculated six accuracy parameters following the procedure in Zhan et al. (2005), and Al-Rawabdeh et al. (2016): user's accuracy (UA; correctness), error of omission, producer's accuracy (PA; completeness), error of commission and overall accuracy (OA) together with the kappa coefficient (k).

Manual landslide event inventory mapping
The manual, visual interpretation of remotely sensed imagery is a common approach in recent landslide inventorying, particularly for mapping landslides produced by one triggering event (e.g. Guzzetti et al. 2004;Malamud et al. 2004;Tsai et al. 2010;Fiorucci et al. 2011;Scaioni et al. 2014) where the landslides are clearly defined. In this study, we produced an inventory of landslides triggered by Hurricane Manuel on September 2013 using satellite images captured in 2011 (7 March 2011) and 2014 (13 April 2014; 16 April 2014 and 12 May 2014) (Image © 2016 DigitalGlobe), i.e. before and after this extreme meteorological event. The visual comparison of these images (i.e. tonal, textural and vegetation characteristics) enabled identification and classification of landslides as debris flows, shallow soil slips and deep-seated features. This method was applied in three main stages: (1) preparation of the preliminary landslide inventory map based on the landslide visual interpretation using satellite images, (2) field work aimed in validation of the preliminary inventory (Guzzetti et al. 2000) and acquisition of data on landslide type and main characteristics in specific areas (Guzzetti & Cardinali 1990), and (3) production of the final landslide inventory map.

Automatic inventory mapping -CCM
We applied a landslide inventorying algorithm, the CCM (Leshchinsky et al. 2015), to supplement manual inventorying and investigate the efficiency of inventorying landslides using automated techniques. Various landslide inventorying algorithms for LIDAR DTMs have utilized image segmentation processes (Van Den Eeckhaut & Herv as 2012), object-based image analysis (Li et al. 2015), the random forest algorithm  or machine learning techniques (e.g. Booth et al. 2009;Li et al. 2015) to automate landslide inventorying. These approaches primarily delineate signatures that are representative of landslide scarps or extents; however, they require adequate training datasets and high resolution for successful mapping. In contrast, CCM detects and uses the shape of topographic features to locate past landslides with minimal topographic parameters representative of landslide morphometry (e.g. headscarps, deposits). This process is performed by first discretizing elevation contours from the LIDAR DTM into evenly spaced nodes to enable calculation of slope between node-node connections on the contours. Using this framework, CCM first identifies landslide scarps based on an upper slope threshold, then finds the steepest connections downslope until a lower slope threshold is reached, mapping the boundaries of a landslide and, ultimately, its extents ( Figure 3). This process demonstrated reasonable accuracy with manual landslide inventories performed using LIDAR DTMs (Leshchinsky et al. 2015), but cited some inaccuracies in correctly identifying headscarps using solely slope as the defining attribute.
To refine the headscarp identification, in this study, the CCM algorithm was modified to allow direct definition of headscarp morphology and subsequent flow analyses to define landslide extents. This process incorporates more flexibility in the landslide inventorying process, enabling input of manually digitized headscarps or automated scarp features to initiate an analysis; however, input of headscarps requires use and combination of site-specific, semi-empirical LIDAR derivatives to highlight regions of perceived headscarp morphology. This process, outlined in Figure 4, requires input of several parameters, including (1) contour intervals, which govern the resolution of the analysis, often selected to be at least three times greater than the raster resolution, (2) on-contour node spacing, also selected to be at least three times greater than the raster resolution, (3) active slope, which governs the minimum slope connections must maintain before terminating, and (4) branching parameter, which governs the number of connections that initiate from a given active node. More details on parameter selection are provided in Leshchinsky et al. (2015).
In the application of CCM to the study regions, the LIDAR DTM was first smoothed with focal statistics by averaging elevations within a 3 £ 3 window to create a smoothed DTM as a preliminary means of filtering noise. Headscarps often exhibit steeper slopes and more notable profile curvature Step 1: Assign contours and nodes.
Step 3: Continue connections from scarp until slope is lower than threshold.
Step 4: Create landslide polygon. than regular hill slopes. To highlight these regions, LIDAR derivatives of slope and profile curvature were obtained through a thresholding technique performed by multiplying slope (S) and profile curvature (C PR ), highlighting perceived headscarp morphology in a Scarp Index, SI, defined as follows: SI values between ¡190 and ¡55 tended to exhibit headscarp features that were deemed appropriate for preliminary delineation of scarp morphology. These rasterized regions were then refined by delineating regions with exposed rock, defined by terrain texture (i.e. roughness representative of exposed rock) as a Boolean to remove headscarp cells. Roughness, R, was defined as the standard deviation of slope using a 3 Â 3-cell moving kernel. This threshold for removing scarps was coincident with regions that have a roughness greater than 1.75 standard deviations, representative of the roughness of exposed bedrock and not headscarp morphology. The remaining scarp pixels were then vectorized by fitting rectangles using the Minimum Bounding Geometry tool in ArcMap Ò . Afterwards, a final refinement was performed by removing scarps with negative average plan curvature. Finally, the Interpolate Shape tool was used to assign unique values to each headscarp, resulting in a set of headscarps that were derived through an automated process.

Logistic regression method
Numerous qualitative and quantitative approaches to map landslide susceptibility exist, employing a wide range of input data-sets to describe landslide-causing factors (e.g. Aleotti & Chowdhury 1999). Qualitative approaches are based on an expert's prior knowledge on the roles of geological and geomorphological factors on landslides. In this approach, an expert develops rules and/or assigns weights to landslide-causing factors applied in mapping susceptibility. Quantitative approaches are based on the integration of landslide-causing factors by using observed statistical relationships with landslides (probabilistic methods), machine learning or by using numerical equations that explain physical mechanism of landsliding (deterministic methods). The utility of LIDAR data in landslide susceptibility mapping has been already proven employing numerous approaches, e.g. logistic regression (e.g. Chen et al. 2013;Jebur et al. 2014;Mahalingam & Olsen 2015;Mahalingam et al. 2016), likelihood ratio (e.g. Chen et al. 2013), frequency ratio (e.g. Pradhan et al. 2012;Mahalingam et al. 2016), weights of evidence method (e.g. Jebur et al. 2014;Mahalingam et al. 2016), discriminant analysis (e.g. Mahalingam et al. 2016), artificial neural network (e.g. Mahalingam et al. 2016), support vector machine (e.g. Jebur et al. 2014;Grejner-Brzezinska et al. 2015;Mora et al. 2015;Mahalingam et al. 2016).
In this study, we followed a probabilistic method of logistic regression proposed by Regmi et al. (2014) using LIDAR-derived landslide-causing factors data-set. The logistic regression approach was selected for two reasons: first, it results in the probability of landslide occurrence in each pixel of the study area, and second, Regmi et al. (2014) showed that the approach is applicable in mapping landslide susceptibility with better accuracy compared to other probabilistic approaches including weightof-evidence and fuzzy logic (Regmi et al. 2010a(Regmi et al. , 2010b(Regmi et al. , 2014. The applied method consists of three main stages ( Figure 5): (1) identification and collection of landslides and landslide-causing factors, together with the conversion of these data into raster GIS formats, (2) landslide classification and sampling, and the extraction of landslide and non-landslide data from maps of landslide-causing factors, and (3) logistic regression employing the data collected in previous stages. The logistic regression uses the following equation to determine the probability of the presence or absence of landslides: where x is the data vector containing multiple independent variables representing the landslide-causing factors, b is a vector of the coefficient(s) of the independent variable(s), and p is the probability of the presence (y = 1) of the dependent variable (landslide occurrence). We performed the calibration of the regression coefficients twice to compare the difference in the factors when derived by manual and automated landslide inventories. During this calibration, multiple factors were initially included and analysed for their contribution. Using both landslide inventories, rainfall-triggered features dominate, which tend to occur along the topographic convergence (Montgomery & Dietrich 1994;Roering et al. 1999;Regmi et al. 2014). Thus, the topographic characteristics are believed to be the most important factors controlling the distribution of landslides in the study area. Variations in geology, land use and soil could also influence the landslide formation; however, for the study area, these data-sets lack the resolution to provide useful landslide susceptibility information and do not show any differentiation. For this reason, the primary landslide-causing factors extracted from LIDAR DTMs were mainly based on topography, hydrology and climate, including elevation, aspect, slope, tangential curvature, plan curvature, profile curvature, distance to stream, flow length, flow accumulation, stream power index (SPI), topographic wetness index (TWI) and solar radiation (see Figures S1 and S2 in Supplementary Material). More details on the selected landslide-causing factors are provided in Regmi et al. (2014).
For each inventoried landslide, we extracted the sample data from the landslide mass centroid. To eliminate bias in the sampling process, we used the same number of landslide and non-landslide cells (the latter were randomly selected across the scene) for the extraction of data from maps of landslide-causing factors. The logistic regression analysis was then conducted using the SPSS Statistics Ò software. The final susceptibility map was produced as a classified probability of the landslide occurrence, which varies from 0 (no susceptibility) to 1 (complete susceptibility). For simplicity, we used four equal intervals of susceptibility: 0-0.25 very low susceptibility (VLS), 0.25-0.5 low susceptibility (LS), 0.5-0.75 high susceptibility (HS), and 0.75-1 very high susceptibility (VHS),  similar to other researchers (Yesilnacar & Topal 2005;Nefeslioglu et al. 2008). The goodness-of-fit for developed models were evaluated based on the significance of the likelihood ratios. Specifically, Cox and Snell R 2 values, Nagelkerke R 2 values and areas under the receiver operating characteristics (ROC) curves (AUC), i.e. the larger the area below the ROC curve, the more accurate the prediction of landslide occurrence (Egan 1975). Using this approach, if all landslides were correctly predicted, the AUC would equal unity.

Manual event landslide inventory
Following Hurricane Manuel, a total of 419 landslides were mapped manually ( Figure 6 and Table 1). These landslides consist of mainly small features, i.e. average length ca. 100 m, width <30 m and area <1800 m 2 (Table 1). Among them were debris flows occurring mostly along zones where topographic convergence predominate (ca. 95% of all identified landslides). These findings are consistent with high values of L/W ratio (up to >10), suggesting fluidization processes produced by heavy precipitation. Observed debris flows are usually long (up to more than 2.5 km) and narrow (32 m in average) features. Larger (average area ca. 18,000 m 2 ) and wider (average width ca. 100 m) deepseated landslides are less frequent (3% of all identified landslides), but did cause significant harm as noted by the deep-seated La Pintada landslide, which covered an area of >120,000 m 2 .

Region No. 1 (North)
We manually inventoried approximately 250 landslides in region No. 1, demonstrating areas between 50 and 130,000 m 2 (Table 1). Surprisingly, a clear correlation between the occurrence of landslides and slope ( Figure 6) was not observed. Generally, landslides are expected to be more common on steeper slopes. Nevertheless, even though important, slope steepness does not seem to play a key role in landslide spatial distribution. For example, the large, deep-seated and catastrophic La Pintada landslide that caused 71 fatalities occurred on a moderately steep slope (25 ) that was lightly vegetated. On the other hand, slopes covered with dense vegetation, regardless of their steepness and orientation, produced mainly narrow and limited debris flows. For example, steep slopes in the densely forested central part of this region exhibited few landslides. The presence of landslide features appears to be related to the distance from rivers and streams (49.2% of landslides were observed in zone up to 200 m from main rivers). In the north-western region, there are quite a few shallow soil slips occurring in proximity to roads. In this region, lithology does not influence the spatial distribution of landslides because the entire study area sits on similar plutonic rocks. No variations in soil development were observed, which is a factor that often influences the spatial distribution of landslides.

Region No. 2 (South)
Landslides identified in region No. 2 are smaller than those is region No. 1 (Table 1). These are mainly narrow, linear features with high values of L/W ratio. The exposure and steepness of slope, vegetation, solar radiation and land use influence their spatial distribution. The number of landslide features on the slope facing south is significantly higher than on the slope facing north, i.e. 142-34, respectively ( Figure 6). The exposure of slope related with the differences in solar radiation determines the use of soil and type of vegetation. Slopes facing south with significantly higher solar radiation are occupied with pastureland and agriculture, thus lacking dense vegetation, making them more prone to mass wasting processes. Alternatively, slopes facing north, with lower solar radiation, are predominantly covered with dense tropical forest that confine and reinforce the natural soils. The differences in manifested landslides with aspect could be also related to steepness. More than 80% of landslides were identified on slopes facing south, which are steeper (mean slope 30.2 ) than slopes facing north (mean slope 24.5 ). The north-facing slopes were mostly devoid of landslides, Figure 6. Landslides identified in both regions using two approaches including manual mapping using satellite images and automatic identification of landslide morphology employing CCM algorithm.
whereas landslides on the south-facing slopes exhibited a tendency to cluster in the steeper sections of these slopes. Similar to region No. 1, the lithology does not influence landslide occurrence because both slopes underlie by the same rock type (i.e. granite to granodiorite), and only the eastern part of the northern slope is made of gneiss.

Automatic landslide inventory
More than 2100 headscarp features were extracted from the LIDAR DTMs. CCM was run using these headscarp features with a contour interval of 3 m, node spacing of 3 m, an active slope of 0.10 and a branching parameter of 2. A total of 1726 landslides were identified, i.e. 1382 landslides (77%) in region No. 1, and 344 (23%) in region No. 2. From this database, we randomly selected 10% of landslides for detailed verification. For both studied regions, we found that ca. 10%-13% of the landslide extents were the result of over-prediction by CCM.

Region No. 1
The CCM algorithm detected 1382 landslides in region No. 1, including the deep-seated landslide in La Pintada (Figure 6). Most landslides were small to medium in size, demonstrating a mean area of landslide extents (measured directly from landslide-extent polygons, not bounding box areas) of approximately 4500 m 2 . Approximately 25% of the identified landslides were below 1000 m 2 . Those smaller landslides occur mainly in the central and northern parts of this region, related with moderate to steep slopes. However, several large landslides exhibited large affected areas, spanning up to 120,000 m 2 ( Figure 6). Those features tended to occur in the southern region, related to the steep topography present there. The identified landslide features primarily occurred in proximity to rivers and streams, but also near roads and other anthropogenic features, such as those in the western and northern parts of the region. The calculated statistics demonstrate a bimodal trend in L/W ratios with a mean ratio of approximately 3.6. Almost 300 landslides demonstrate L/W ratios less than two, typical to slumps and deep-seated landslides that do not fluidize and run out significantly. However, many of the remaining landslides demonstrated large L/W ratios (up to >10), suggestive of fluidization, consistent with debris flows and mudslides that occur after heavy precipitation.  Figure 6). We observed a significant difference in the number of landslide features based on the orientation of the slopes. On the north-facing slopes, only 73 landslides were identified (21.2% of all landslides in this region), whereas on the opposite slope facing south, we recognized 271 features (78.8%). Landslides observed in this region are generally smaller than those in region No. 1, realizing average landslide areas of less than 70% of those in region No. 1 (Table 1). L/W ratio shows similarly high values as in region No. 1, suggesting predominance of debris flows and mudslides that occur after heavy precipitation.

Accuracy assessment
We assessed the accuracy of both landslide inventories by comparing the obtained results with the reference data validated in the field ( Table 2). The results indicate that landslide inventories produced by both the manual interpretation of satellite images and automatic CCM algorithm show very high values of overall accuracy, i.e. >90% (99.67% and 90.65%, respectively). The producer's accuracies (completeness) and user's accuracies (correctness) were found to be 99.04% and 99.08%, respectively, for the manual approach, and 75.74% and 75.47%, respectively, for the automatic approach. Correctness and completeness of manual inventory are also evidenced by very small values of errors of commission and omission (<1%; Table 2). Small incorrectness is probably related only with errors in the visual and manual digitizing process related with horizontal inaccuracy present in the satellite images, as well as potential vegetation overgrown that could cover some parts of shallow landslides. Higher error values calculated for the CCM approach (Table 2), likely related to over-prediction of landslides related to the primary objective of this method, which is to capture all past landslides. Over-prediction is also particularly evident in river valley bottoms and anthropogenic forms where the CCM approach identified some large rocks, Cohen's Kappa coefficient ranging from 0 to 1 provides a measure of the degree of agreement of two landslide inventories (Cohen 1960). Kappa value varies between two landslide inventories, from 0.51 for the CCM automatic inventory to 0.98 for the manual identification of landslide features (Table 2), indicating good and strong agreement between the extracted results and the reference data, respectively.

Logistic regression
We produced maps of landslide susceptibility using the logistic regression method based on landslide-causing factors data derived from the LIDAR DTMs and the two landslide inventories. The goodness of fit of the developed models is presented in Table 3.  Figure 7). For example, the model based on event inventory shows high and very high susceptibility zones, defined as probability of landslide occurrence >0.5, mainly related to regions of high topographic convergence (Figures 7 and 8(a-c)). This agrees with the predominance of debris and mudflows triggered by extreme precipitation related to Hurricane Manuel. On the other hand, the model based on CCM inventory shows high probability and susceptibility on the steep sections of southfacing slopes (Figures 7 and 8(b,d)). In some regions, the model based on the manually derived event inventory overestimates the landslide-affected areas as evidenced by high and very high landside susceptibility predicted for the area of La Pintada village (Figure 8(a)). Similar over-prediction was not observed for the model based on the CCM inventory (Figure 8(b)). Apart from differences described above, both models present higher probability and susceptibility to landslides for the southern part of the region No. 1 than for the central and northern regions (Figure 7). Susceptibility on slopes facing west and east do not show a notable difference; however, higher susceptibility is observed for east-facing slopes. The areas under the ROC curve indicates similar overall prediction accuracy for both models, i.e. 0.77 for model based on the manual event landslide inventory and 0.78 for model based on the CCM inventory (Figure 7 and Table 3).
The spatial distribution of susceptibility to landslides is strongly influenced by the following landslide-causing factors: TWI, proximity to stream, slope steepness, SPI and tangential curvature (Table 4). Additional factors, such as elevation and aspect, were also important for the susceptibility model based on the CCM inventory. TWI is an effective parameter in analysis of susceptibility to landslides triggered by heavy rainfall, because it enables determination of areas that are most likely to become saturated during heavy precipitation (Moore et al. 1991). Thus, high coefficient of TWI is considered reasonable as this region exhibited a propensity towards debris flows triggered by exceptionally heavy precipitation (e.g. as demonstrated by the aftermath of Hurricane Manuel). Similar dependence was not observed for the susceptibility model based on CCM inventory, as it includes also other landslide types. SPI measures the erosive power of water flow (Moore et al. 1991), where higher values are usually found at the foot of a slope. The high influence of this factor together with the distance to stream on the produced susceptibility models indicates that the landslides are more probable to occur closer to the streams, at the foot of the slope.

Region No. 2
In region No. 2, there is a pronounced difference between landslide susceptibility based on slope orientation, regardless the type of landslide inventory used to calibrate the logistic regression coefficients (Figure 9). Steep slopes facing south present significantly higher landslide probability than north-facing slopes. Furthermore, areas of high susceptibility cluster primarily on the south-facing slope.
The areas under the ROC curve indicate 0.87 and 0.83 overall prediction accuracy (AUC) for susceptibility models based on the manually derived event inventory and CCM inventory, respectively ( Figure 9 and Table 3). As shown in Table 4, the most important landslide-causing factors used in Table 3. Statistics of the logistic regression applied to produce maps of landslide probability and susceptibility using two different inventories.
Likelihood ratio test the logistic regression analysis are hydrological parameters such as TWI and distance to stream, and hillslope morphological parametersslope steepness, but also profile curvature for the CCM-based model ( Table 4). The strong impact of TWI is related to the predominance of debris flows and mudslides that occur after heavy precipitation. A significant influence of distance to stream on the produced maps of susceptibility suggests that slopes undercut by river incision tend to be an important triggering factor for landsliding processes, as observed in many regions (e.g. Regmi et al. 2014). During Hurricane Manuel, some river discharges in Sierra Madre del Sur was up to 900 times higher than average values from 1955 to 2012, exhibiting stream levels higher than 7 m (CONAGUA 2016). These exceptionally high flows and water levels likely produced many landslide features. Again, the steepness of slope in the study region exhibits significance in landslide activity, demonstrated by the larger number of landslides on steep, south-facing slopes versus gentler north-facing slopes.

Validity test
We tested the validity of presented landslide susceptibility analyses using the approach by Regmi et al. (2014). First, from both inventories we randomly selected training and testing data-sets (both with 50% of landslides) using the Subset Features tool in ArcMap 10.2.2 Ò (http://desktop.arcgis.com/en/). Consequently, we produced maps of landslide susceptibility based on the training data-set. We then verified the validity of these models by comparing the probability to landslides with landslides from the testing database. The resulting high AUC values for both study regions and both inventories (0.75-0.82; Figure 10) suggest that the developed models present reasonable accuracy. We also validated the susceptibility maps based on one landslide inventory (either landslide event inventory or automated CCM inventory) by comparing the results with spatial distribution of landslides from the second inventory. This was done by extracting values of landslide probability to points representing landslide mass centroids, and consequently producing ROC curves. Obtained values are higher for the region No. 2 than for the region No. 1, which agrees with overall prediction accuracy of developed susceptibility models (Figures 7, 9 and 10, and Table 3). Susceptibility models produced based on CCM inventory predicted landslides triggered by Hurricane Manuel with high accuracy (0.70 for the region No. 1 and 0.86 for the region No. 2; Figure 10). On the other hand, models produced based on the event inventory predicted with high accuracy the older landslides identified by the CCM algorithm (0.77 and 0.81 for regions No. 1 and No. 2, respectively, Figure 10).

Comparison of landslide inventories
In comparing the inventories (Figure 6), CCM identified approximately four times more landslides than the visual interpretation of satellite images. Schulz (2007) obtained a similar difference between the number of landslides identified using LIDAR data and aerial photographs, related mainly to the differences in data resolution between the two techniques. In this study, the relatively high discrepancy in the absolute number of landslides and some landslide characteristics is strictly related to the differences in the objectives of the inventorying procedures. Whereas visual interpretation of satellite images was aimed to detect landslides triggered by the Hurricane Manuel, the CCM algorithm was employed to identify all past landslide features by analysing the morphology of the area represented by the LIDAR DTMs. Thus, CCM captured landslides produced by this hurricane, but also older features, such as those produced by rainfall related to other extreme meteorological events or by other factors, such as earthquakes or anthropogenic influences. High accuracy values, comparable or even higher than those reported in the literature prove the validity of both inventories (see Table 2). On the other hand, comparison of two inventories, performed by superimposing manually mapped landslides over CCM features revealed that, in some cases, the modelled landslide extents do not capture the entire run-out extents present on the manual inventory map. This is likely due to the shallow gradients of the landslide repairs at its toe, preemptively terminating the flow of the CCM algorithm ( Figure 6). Even though the results obtained by both methods show a discrepancy in the number of identified landslides (Table 1), the general pattern of landslide spatial distribution is similar (Figure 6). Both inventories show that most landslides in the region are small, but some large slope failures exist. In general, landslide features identified by CCM are two times longer and two times wider than those mapped manually (Table 1). Landslides are distributed randomly in the region No. 1, whereas their distribution in the region No. 2 shows a clear difference dependent on slope aspect, i.e. four times more landslides were identified on the south-facing slope than on the opposite one. This observation is perhaps related to solar radiation, slope steepness and respective vegetation density. The central and southern parts of region No. 1 show the most significant discrepancy between landslide inventories. Event inventory shows that extreme precipitation related to Hurricane Manuel did not produce large landslide features in this region, but mainly narrow debris flows. However, it is likely that large failures could have happened in the past during more intense storms than Manuel. This is evidenced by the existence of large number of large-sized landslides identified by CCM that are not present in the manual inventory.
The landslide inventory produced by CCM demonstrates more exactness, primarily because it identified features four times smaller in area than the smallest landslide mapped based on satellite images (Table 1). This could be related to the large resolution differences between the data sources used in the manual and automatic landslide identification.
Landslides triggered by the heavy precipitation that occurred during Hurricane Manuel were primarily debris flows and mudslides, i.e. fluidized features that demonstrate long run-out lengths and rapid movements. Similarly, many of the landslides identified by CCM show highly elongated shape (high L/W values) typical to slides that travel long and very long distances. This suggests that most of the landslides identified by CCM could have been also triggered by exceptional precipitation related to extreme meteorological events. Hence, it is likely that rainfall-triggered landslides predominate in studied regions. This finding is of particular hazard concern as fluidized landslides and debris flows move quickly with little warning.

Susceptibility to landslides and landslide inventories
The applied logistic regression method and high-resolution LIDAR data enabled production of high-resolution landslide susceptibility maps in the Guerrero mountains, with reasonable prediction accuracy, ranging from 0.77 to 0.87, depending on the region and type of landslide inventory used for the calibration of regression coefficients (Figures 7 and 9, and Table 3). Similar AUC values obtained for models based on manual and automated landslide inventories suggest that when using high-resolution LIDAR data, the landslide susceptibility mapping process can be automated without a loss of accuracy. Furthermore, the susceptibility model based on automated inventorying might demonstrate even higher prediction accuracy than that based on manual landslide inventorying (compare AUC values for region No. 1 in Table 3).
Most of the previous studies using logistic regression method for susceptibility mapping show accuracy of resulting models ranging from approximately 0.8 to 0.9 (Table 5), of which the highest values were obtained for small areas with only few landslides (<20) (Nefeslioglu et al. 2008). Data presented in Table 5 suggest that the fewer landslides employed in the logistic regression calculations, the higher the prediction accuracy of the resulting model. However, susceptibility models based on fewer landslides might result in high, but actually false non-significant prediction accuracy.
In this study, approximately 2000 landslides (CCM inventory) and >400 features (manual event inventory) were used for the calibration of regression coefficients (Table 1), exhibiting reasonable accuracy. Similar values were obtained by Regmi et al. (2014), who also used large landslide database and followed the same logistic regression analysis procedure. Hence, reasonable accuracy can also be attained using larger inventories, possibly expanding applicability to larger regions. Moreover, calculations were performed from both inventories, regardless of their type or size, whereas previous studies realized higher prediction accuracy for only one type of landslide produced by a singular triggering event (e.g. Regmi et al. 2014). This suggests that already high AUC values could be improved.
Previous studies using logistic regression approach for landslide susceptibility maps, with a few exceptions (e.g. Chen et al. 2013;Olsen et al. 2015) based mainly on DEMs with horizontal resolutions greater than 10 m (Table 5). Here, LIDAR DTMs with 1 m horizontal resolution were used, significantly improving the resulting susceptibility models. Even though the obtained values of prediction accuracy (AUC) are comparable with data in the literature (Table 5), the models based on LIDAR data enabled enhanced precision in the identification process of zones with high and very high susceptibility (Figures 7-9). The resulting models show clear differentiation of landslide probability and susceptibility along and across slopes of detailed scales based on the variations of topographic, hydrologic and climatic landslide-causing factors. Attempts to map landslide susceptibility using LIDAR data already have demonstrated applicability, but have generally alternative approaches than logistic regression models (e.g. Chigira et al. 2004;Schulz 2007;Jaboyedoff, Choffet, et al. 2012;Jaboyedoff, Oppikofer, 2012;Pradhan et al. 2012;Chen et al. 2013;Jebur et al. 2014;Grejner-Brzezinska et al. 2015;Mahalingam & Olsen 2015;Mora et al. 2015;Mahalingam et al. 2016). In this study, susceptibility models with reasonable accuracy were attained, even though causative factors for landsliding were only extracted directly from the LIDAR DTMs. This indicates the utility of high-resolution LIDAR data for performing susceptibility analyses.
Susceptibility models based on different landslide inventories can lead to contrasting patterns in spatial distribution of landslide probability and susceptibility, despite the fact that both show similar AUC values (Figures 7 and 8). Significant differences in susceptibility models are related to the different goals of each inventorying procedure, and to the type of landslides used in the logistic regression analysis. The event inventory captured mainly debris flows and mudflows triggered by extreme precipitation related to Hurricane Manuel. Those landslides occur mainly in zones of topographic convergence (Montgomery & Dietrich 1994, Roering et al. 1999Regmi et al. 2014). For this reason, the susceptibility model based on this inventory shows high and very high landslide probability and susceptibility in regions of topographic convergence, i.e. gullies, deep stream valleys, etc. (Figures 7  and 8). This might lead to over-prediction as evidenced by high and very high susceptibility zones identified by the logistic regression approach in the La Pintada village (Figure 8(a)). On the other hand, model based on CCM inventory shows high probability and susceptibility for all types of landslides. The CCM inventory, in contrast to the event inventory, identified, apart from debris flows and mudflows, a large number of deep-seated features ( Figure 6). Their presence and spatial distribution strongly influenced the susceptibility model, which does not focus only on zones of topographic convergence but also focuses on the slope curvature, aspect and distance to streams (Table 4 and Figures 7 and 8(b,d)). Thus, the resulting susceptibility model based on CCM inventory demonstrates broader applicability, showing the general probability and susceptibility to landslides regardless of their type and triggering mechanism, whereas the model based on event inventory presents the probability of occurrence of landslides triggered only by one event, e.g. extreme precipitation. This concept has been demonstrated in other cases when applying event inventories in susceptibility modelling that predicted mainly landslides of the same type (e.g. Lee et al. 2008;Deb & El-Kadi 2009). In summary, our results demonstrate that the entire process of landslide susceptibility mapping can be automated by combining CCM landslide inventory and logistic regression approach using high-resolution LIDAR data. However, further studies are needed to improve landslide susceptibility models by applying additional or different landslide-causing factors, using different sampling techniques, and landslide classifications.
Good agreement is achieved when comparing the inventories with susceptibility maps, with many landslide features observed in areas of high and very high susceptibility (landslide probability >0.5). This is valid for both landslides triggered by Hurricane Manuel and landslide extents identified by CCM ( Figure 10). Locations of features identified by CCM agree with zones of high and very high susceptibility based on susceptibility modelling using landslides triggered by Hurricane Manuel. Similarly, landslides included in the manual inventory occur mainly in the high and very high susceptibility zones determined using the CCM inventory. This indicates that the susceptibility models are effective and compounds the predominance of rainfall-induced landslides within the studied regions. This finding is corroborated by significance of landslide-causing factors related to topography (i.e. slope steepness and tangential curvature) and rainfall (i.e. TWI, SPI and proximity to stream) in the presented models of landslide susceptibility.

Conclusions
We tested and compared two different methods for landslide inventorying and evaluated their effect on landslide susceptibility mapping. The automated CCM enabled a fourfold increase in inventoried landslides in comparison to manual interpretation of satellite images, increasing exactness of the analysis with the use of high-resolution LIDAR DTMs. This suggests that CCM may present a useful supplement to manual inventorying, particularly due expedience and objectivity implicit in an automated analysis; however, LIDAR topographic data are required for implementation. The results from CCM demonstrate good performance for inventorying large to medium landslides, but demonstrate difficulty in finding very small slumps due its discretization procedure. Finally, the CCM inventoried most of the landslides that occurred following Hurricane Manuel as well as many older features, exhibiting utility in not only finding recent slope failures, but also older, weathered features, producing complete landslide inventory.
We demonstrated that susceptibility models based on an automated landslide inventory show high overall prediction accuracy (up to 0.83) comparable with models based on manual event landslide inventory, contingent upon the availability of high-resolution LIDAR data. Thus, landslide susceptibility hazard maps can be made automatically using the proposed-here techniques: (1) CCM algorithm to identify landslides based on the shape of topographic features, and (2) the logistic regression method to produce map of probability and susceptibility to landslides. The results of the combination of both methods were verified and compared with actual existing landslides mapped using manual interpretation of satellite images.
Landslide susceptibility modelling was performed using a logistic regression approach, demonstrating applicability for the tropical mountains of Guerrero in south-west Mexico. We achieved satisfactory accuracy, demonstrated by AUC values between 0.77 and 0.87 for overall prediction accuracy, depending on the region and inventory used for the calibration of regression coefficients. The resulting models identified the majority of mapped landslides (both manual and automated inventories) within the high and very high susceptibility zone (defined as probability of landslide occurrence >0.5).
Automated remote sensing techniques for mapping landslides and susceptibility to landslides demonstrate promise for expanding the spatial scales of planning and disaster preparedness efforts. However, these approaches are constrained by the availability and high cost of high-resolution LIDAR data, in particular in developing countries, required to reach sufficient accuracy of the resulting products. These findings provide further justification for making LIDAR and other high-resolution remote sensing data more accessible through freely access to available data and sharing systems, such as Open Topography (http://www.opentopography.org/) initiative. Fortunately, the increasing utility of LIDAR technology with advances in computational capabilities demonstrates promise for future data collection and economic feasibilitya promising tool for landslide preparedness.
Special acknowledge to Mario Paredes for his contribution to this project and Nathan Lyons for the initial landslide inventory and data collection. We thank the editor and the anonymous reviewers for providing comments that helped to improve our manuscript.

Disclosure statement
No potential conflict of interest was reported by the authors.