Different Approaches to Use Morphometric Attributes in Landslide Susceptibility Mapping Based on Meso-Scale Spatial Units: A Case Study in Rio de Janeiro (Brazil)

Landslide susceptibility maps are widely used in landslide hazard management. Although many models have been proposed, mapping unit definition is a matter that still needs to be fully examined. In the literature, the most reported mapping units are pixels and slope units, while in this work, developed in the Rio de Janeiro region (Brazil), the use of drainage basins as a mapping unit is examined; even if their use leads to the definition of maps with a coarser spatial resolution than pixels-based maps, they convey information that can be easily and rapidly handled by civil defense organizations. However, for the morphometrical characterization of entire basins, a standardized procedure does not exist, and the susceptibility results may be sensitive to the approach used. To investigate this issue, a random forest model was used to assess landslide susceptibility, using 12 independent variables: four categorical (land use, soil type, lithology and slope orientation) and eight numerical variables (slope gradient, elevation, slope curvature, profile curvature, planar curvature, flow accumulation, topographic wetness index, stream power index). For each basin, the numerical variables were aggregated according to different approaches, which, in turn, were used to set up four different model configurations: i) maximum values, ii) mean values, iii) standard deviation values, iv) joint use of all the above. The resulting maps showed noticeable differences and a quantitative validation procedure showed that the best configurations were the ones based on mean values of independent variables, and the one based on the combination of all the values of the numerical variables. The main outcomes of this work consist of a landslide susceptibility map of the study area, to be used in operational procedures of risk management and in some insights on the best approaches to aggregate raster cell data into wider spatial units.


Introduction
Landslides are a geomorphological process responsible for deaths and damages worldwide [1] and landslide susceptibility maps (LSM), which represent the spatial probability of landslide occurrence in a given location [2], are one of most widely used tools for the related hazard management. Landslide susceptibility is usually assessed by a study of the spatial distribution of static predisposing factors and their relationship with the location of the landslides occurred in the past [3][4][5], as it is postulated that future slope failures will be more likely to occur under the conditions that led to past and present instability [6,7].
LSM have been produced applying a wide range of mathematical techniques, from the most traditional statistic approaches like frequency ratio [8], discriminant analysis [9,10] and logistic regression [11][12][13], to more recent and more advanced techniques, like artificial neural network [5,14], machine learning [15,16] and multi criteria decision analysis [17]. Many studies highlighted that different models can produce different results but, to date, no model has been proved to be better than all the others [5,8,10,12,16,17]. Indeed, another series of studies showed that the same model in the same test site can provide very different outcomes, depending on the approaches used to prepare the input data [18][19][20].
Independently from the mathematical model used in the landslide susceptibility assessment, geographic information systems (GIS) remain a fundamental tool for LSM studies, as they can be used in data preparation, data analysis and in drawing the outputs [21][22][23]. One of the advantages of the use of GIS in LSM is the possibility of overlaying different thematic maps and analyzing the spatial distribution of their attributes across appropriate spatial units. Concerning the spatial units of analysis and of representation of the resulting outputs, a recent review by Reichenbach et al. [24] showed that pixels are by far the most used mapping unit in the international literature, and only a minor number of works exists that use unique condition units (units defined in terms of the possible combinations of a limited set of input parameters in a given space), slope units or other techniques of segmentation of the territory [13,[25][26][27][28].
As a consequence, we believe that the application of unit-based landslide susceptibility models instead of pixel-based models has not been deeply investigated. For instance, the recent literature is rich of examples of comparisons among different models [8], sensitivity analyses to model settings [19], sampling strategies [29], parameterization of variables [30] or pixel cell size [31], but despite a few attempts [32,33], scarce attention has been paid on the different approaches to generalize input data from the original raster format to the polygonal nature that characterizes other spatial units.
The main objectives of the present work are twofold. The first practical objective is to produce an LSM for the mountain region north of Rio de Janeiro, Brazil, to assist civil defense agencies in hazard management procedures. Another important research objective of this work is to cope with a source of uncertainty that affects LSM based on spatial units larger than the pixel size of the input parameters: how to generalize the fine-level pixel information at the level of the chosen spatial unit. In this study, a state-of-the-art machine learning algorithm (random forest [19,34]) has been used for landslide susceptibility modeling, and the hydrologic basins of the State of Rio Janeiro have been used as spatial units for the analyses. The approaches of characterizing the basins in term of the mean, the maximum and the standard deviation of the values of the numerical predisposing factor has been tested, highlighting noticeable differences in the resulting LSM and providing insights to decide which is the best option and if these strategies can be combined.

Test Site Description
The study area comprises the municipalities of Petrópolis, Teresópolis and Nova Friburgo (Figure 1), which, together, occupy an area of 2500 km 2 . The climate of the region is classified as Cwb (altitude subtropical climate), according to the Köppen classification, resulting in dry winters and rainy and fresh summers. The study area is located in the Serra do Mar complex, an important mountain chain that extends for approximately 1500 km along the southeastern coast of Brazil. The geology of the area is associated with a wide range of folds (Proterozoic age), composed predominantly of high-grade metamorphic rocks (gneisses) with well-defined foliation in the SW-NE direction and fractures in several orientations. Associations of syn-tectonic igneous rocks (granitoids) are also common [35]. Petrópolis and Teresópolis are located in a particular Serra da Mar portion, called the Serra dos Órgãos, which consists of a mountainous wall, raised by tectonics, which delimits the bottom of the Guanabara Bay basin. The unit is characterized by an abrupt orographic barrier with WSW-ENE direction. The highest portion of the Serra dos Órgãos, between Petrópolis and Teresópolis, is a plateau raised to over 2000 m of altitude, which contrasts with the extremely rugged relief of the surrounding mountain escarpments, characterized by very steep slopes and sharp tops. Between Teresópolis and Nova Friburgo, the Serra dos Órgãos maintains its imposing monolithic aspect, with altitudes between 1100 and 1300 m to the west, and 1400 to 2000 m to the east.
At lower altitudes, the relief is smoother, with intramontane hills characterized by convex-concave slopes and rounded tops, with colluvium and alluvium sedimentation at the footslope. In these areas, the regolith is composed by thick saprolitic and colluvial deposits that together can occasionally reach up to 50 m of thickness. The valley bottoms are narrow and develop along persistent tectonic fractures, in which only the larger rivers generate fluvial deposits. Near the valleys, there are escarpments with rocky outcrops and steep slopes, higher than 35 degrees. Talus and colluvium deposits with abundant rock blocks are usually present at the base [36].
The original vegetation was tropical forest, named "Mata Atlântica", but the human occupation and the need for areas available for planting and animal husbandry substantially reduced the natural vegetation. Therefore, areas with high slope, which, at first, were occupied by primary forests, became unstable after the removal of natural vegetation [37]. Also, in the 1960s, the urban area began to expand, occupying places with steep slopes, and consequently, houses were built in hazardous areas [38].
Mass movements of the study area are characterized mainly by rapid shallow translational landslides and debris flows [37]. Typically, the depth of the landslides is between 0.5 and 2.0 m, within the saprolites, in the slope colluvium or at the soil-rock interface [37]. In the study area, shallow translational slides and debris flows are triggered by rainfall. In some cases, they are combined together (e.g., complex landslides formed by translational slides evolving into debris flows) and are influenced by very similar controlling factors; indeed, in the landslide susceptibility literature, these kinds of phenomena are frequently analyzed together [10,39].
In January 2011, 11th and 12th, this region experienced one of Brazil's biggest disasters, as it rained about 320 mm in 48 hours, leading to hundreds of landslides [40]. This event caused 918 deaths (officially) and the 3 studied municipalities were the most affected. Currently, these municipalities are monitored by the National Center for Monitoring and Early Warning of Natural Disasters-Brazil (CEMADEN) 24 hours per day, and receive warnings for landslides and flooding. It is therefore necessary to constantly update the mappings of the susceptible areas, in order to improve the alerts of natural disasters, and to direct attention to areas with the highest potential for landslides.

Random Forest
To generate the landslide susceptibility maps, in this work, the random forest model has been used. The random forest (RF) model is a nonparametric, multivariate, machine learning technique, which was proposed by Breiman [34,41] and first used in landslide susceptibility studies by Brenning [42]. Since then, it has rapidly gained widespread consolidation through many researches and case studies, as it is considered a relatively powerful approach in classification, regression and unsupervised learning [16,[43][44][45]. RF contains a series of binary tree predictors, which are generated by using a random selection of the explanatory variables (or independent variables) to split each binary node (yes/no), and to perform a classification of the target dependent variable (in LSM studies, the presence or absence of landslides). Some of the observations are used for internal testing to evaluate the predictive capability of each predictor tree. This information is used to iterate the procedure hundreds of times by growing other random trees (hence the name "random forest"), and to iteratively adjust the prediction effectiveness. Once the best predictor tree is identified, it is applied to the whole study area, to define the LSM. RF has a great predictive performance, and runs fast by summarizing a large number of classification trees, especially when dealing with large amounts of data.
The random forest algorithm used in this work includes an internal test procedure, which can be used to assess the relative error that would be committed by the model, in case an independent variable is not used (out of bag error (OOBE)). This feature allows for ranking the variables according to their importance, and to implement a procedure of forward selection of parameters, in which the redundant, pejorative or uninfluential parameters can be automatically identified and discarded. Comparing to methods with single classifier combinations, the random forest has several advantages: both categorical and numerical variables can be used simultaneously; the algorithm is deemed suitable for analyzing nonlinear variables without considering the matter of multi-collinearity; it is quite robust with respect to overfitting problems; no particular statistical distribution of the data is required [19,34].

Model Configuration
The study area was divided into drainage basins, using the digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM), modified by Valeriano [46], with 30 m spatial resolution, and the ArcHydro tool, present in the ArcGIS software, totaling 1253 basins, with an average area of 1.9 km 2 . The choice of this spatial unit is mainly linked to the operational procedures in which the resulting map is planned to be used. As stated before, pixel is by far the most used spatial unit for landslide susceptibility analysis and mapping [24], and it has the advantage of a very fine spatial resolution. However, the pixel output is very complicated to understand and to exploit for untrained personnel, as it carries a very fine spatial detail that is difficult to properly address in the operative response. Basins are coarser spatial units, but they have the advantage of referring to a portion of the territory that is easy to understand and in which a unitary counteraction can be implemented in case of alarm.
The random forest model used as dependent variable the landslides and as independent variables a series of physical parameters (land use, soil type, lithology and slope aspect as categorical variables and slope gradient, elevation, slope curvature, profile curvature, planar curvature, flow accumulation, topographic wetness index, stream power index as numerical variables), as reported in Table 1. For the landslide data, the information used for the model was the presence or absence of landslides (binary classification between basins affected by landslides or not). Regarding the independent variables, different approaches were pursued to upscale the information to the basin scale. For categorical parameters, the most frequent class of each basin was considered as the most representative (majority approach); for numerical parameters, different variables were created considering the maximum value, the mean value and the standard deviation of the values encountered in each basin. For instance, the physical parameter "slope gradient" was used to derive three possible independent variables, namely "mean slope" (the mean value characterizing each basin), "maximum slope" (the highest value encountered in each basin) and "standard deviation of slope" (the standard deviation of the values encountered in each basin). A crucial step in LSM analysis is the approach used to sample the variables to train and validate the model. As in any other statistical procedure, the size of the dataset influences the results, therefore the higher the number of samples to perform the statistical calibration/validation of the model, the more reliable the results obtained [19]. The same authors [19] demonstrated that a random sampling improves the predictive capability of the map and that the susceptibility model should also be trained/validated with respect to information about non-landslide locations, to avoid a generalized hazard overestimation. Regarding the proportion between the calibration and dataset samples, it is common practice to split them according to a 20/80 or 30/70 ratio [5,8,25,26,37,38]. As a consequence, in our application, all landslide basins were included in the analysis dataset. Afterwards, to get an unbiased susceptibility assessment, SPSS software was used to randomly select for the analysis also a number of non-landslide basins equal to that of landslide-affected basins. From the whole set of used basins (50% with landslides and 50% without landslides), 432 (70% of the dataset) basins were randomly selected (with SPSS software) for model training, and 196 (30% of the dataset) basins were selected for validation.
To investigate the research questions posed in the introduction, four susceptibility assessments were performed separately, using four different configurations of input parameters. The difference in these configurations consists in the method to consider numerical variables in each unit basin: (i) MEAN VALUE, (ii) MAXIMUM VALUE, (iii) STANDARD DEVIATION, (iv) ALL the three above were used conjunctly. The motivations behind this series of tests is that all these methods have some points of strength: mean values would provide the best possible "description" of the average morphometric characteristics of each basin; maximum values would highlight the most extreme morphometric situations, which are often more linked to the initiation of landslides than averaged representations; standard deviation of morphometric attributes could be used as indicators of the variability of the landforms. The workflow of the research is summarized in Figure 2.

Landslide Inventory
The availability of a complete and accurate landslide inventory is a fundamental requisite for landslide susceptibility assessments [47]. Traditionally, landslide inventories can be compiled using different sources of information including remote sensing, field surveys, internet, official report events [48][49][50][51]. There is no single database of natural disaster events in Brazil, so different data sources have been used, and the collected data have been organized and homogenized.
In this study, landslide data from 2008 to 2018 were obtained from scientific papers, event reports provided by CEMADEN, media, or mapped using GIS and satellite images. Overall, 3483 landslides were inventoried in a GIS-based geodatabase used for model calibration and validation: 2930 in Nova Friburgo, 348 in Petrópolis and 205 in Teresópolis. Data of the January 2011 landslide event for Nova Friburgo were provided by the municipal civil defense, which were mapped using satellite images available by the Google Earth in 01/19/2011 and visual interpretation ( Figure 3). The scars were mapped along their entire length with polygon geometry. Landslide occurrences during the same event in Petrópolis and Teresópolis have been mapped by visual inspection of the Google Earth images acquired in 04/11/2011 and 04/18/2011, and were mapped as points positioned where the slide started. Other events that occurred before and after January 2011 were obtained from CEMADEN, which uses data provided by the local civil defense and media, and also searches for landslides occurrences on the internet. In the aftermath of an event, the civil defense provides CEMADEN with a report containing the localization of the landslides. This information was used to include each landslide in the database as a point with a varying degree of precision (exact coordinates or approximate location). For Nova Friburgo, landslide occurrences found in the work of Oliveira [52] were also used, which were provided by the Nova Friburgo civil defense from November 2008 to December 2011. These different sources provide very heterogenic information about landslides: sometimes they are mapped as polygons, sometimes as points. In addition, polygons may represent the only depletion zone or the whole landslide (from depletion to accumulation zone), while points sometimes are referred to the triggering area, and sometimes to the impact with infrastructures. Lastly, the spatial accuracy of the mapped elements is also not homogeneous: many landslides are mapped with high spatial accuracy, but in some occurrences, the spatial approximation is in the order of some meters. However, this level of uncertainty and inhomogeneity is not a problem if drainage basins are used as a spatial unit for analysis and for mapping, since, for each landslide, the spatial information needed is the basin in which it occurred, and there is no need to exactly pinpoint its location as in pixel-based LSM analyses. All the landslides were integrated in a geo-database in a GIS system, and drainage basins were classified as landslide-basins and no-landside-basins. To assess whether a basin is affected by landslides or not, the spatial detail of the elements of the geodatabase is sufficiently accurate and the mapping approach (whether the scars or the landslide bodies were mapped by polygons or points) is uninfluential, as, usually, each landslide evolves within the same basin.

Explanatory Variables
The physical parameters and the derived explanatory variables (aggregated at the basin scale) used in the susceptibility model are listed in Table 1. The aggregation technique varies, depending on the type of the variable. For the categorical parameters, the majority class present in the entire basin was considered. For the numerical parameters, the mean, the maximum value and the standard deviation were calculated on each basin, thus generating three different variables from each parameter. A brief overview of all the parameters is reported hereafter.
The land use and land cover (LULC) map is a categorical parameter, and was created using Landsat TM5 images, which have 30 m spatial resolution. Four images were used to cover the entire study area: orbit 216, scenes 075 and 076 (09/04/2011) and orbit 217, scenes 075 and 076 (08/28/2011), obtained by the National Institute for Space Research (INPE), are shown in Figure 4. To define the LULC map, the year 2011 was selected, because most of the landslides occurred in this year. The method chosen to make the classification of land use was the neural network tool provided by the ENVI software. This technique replicates the neural structure of the human brain, and can be trained to recognize classification patterns in a complex scenario, such as a satellite image. Pixel samples were selected for the urban area, pasture, eucalyptus, secondary forest, primary forest and water classes ( Figure 5). The omission and commission errors in the final map was corrected using raster editing and high-resolution spatial images available in Google Earth™. Land cover is important for studies of mass movements in the sense that the type or absence of vegetation cover directly influences the soil infiltration capacity [35,53,54]. It is evident that the degradation of vegetation cover coincides with widespread landslides of the slopes, but there are also cases of large landslides trigged by heavy rains in regions covered by forests. As the natural vegetation has been removed and replaced by other land covers, the weakening and destabilization of the slopes occurs [55]. Three different types of forest were mapped, because they have a very different structure, with different implications in slope stability. The primary forest is denser, and the root system is more compact and can better sustain the soil. The secondary forest is sparser, and is the forest in regeneration stage. The roots system is not well-formed and does not support the soil in the landslide events very well. The eucalyptus forest in Brazil is a non-native forest, and is planted for commercial purposes (cellulose, wood and other products). The large spacing between the trees can be a disadvantage for hillslope resistance to landslides.
The soil type parameterization was derived from the Soil Map provided by the Rio de Janeiro Environmental Secretariat, at the scale 1:100,000 ( Figure A1). Different soil types have different physicochemical characteristics and, consequently, varying degrees of susceptibility in relation to landslides. The soil types more connected to shallow landslides are little cohesive, shallow, colluvial and porous [56]. Unstructured soils tend to be more common in steeper slopes or in talus deposits.
Lithological classes were derived from a lithology map provided by the Brazilian Geological Survey (CPRM), at the scale 1:100,000 ( Figure A2). The most common lithology types present in the study area are shown in Table 1. Lithology is responsible for the composition, texture and degree of resistance of rocks and soils, as well as other factors, such as permeability and shear strength, which directly influence the susceptibility to landslides [30,57]. Different rocks have different degrees of stability. All the other parameters were derived from the digital elevation model (DEM) in a geographic information system software. The DEM was obtained by CPRM database and was generated from SRTM or Aster data, with spatial resolution of 30 m, and represents the elevation map ( Figure 1). The aspect ( Figure A3) indicates the slope orientation, and it was selected because in the study area the slopes facing south and southeast are those receiving more rain during the arrival of a cold front from the south, which is the condition that causes more landslides in the country.
The slope angle ( Figure A4) is probably the most used parameter in landslide susceptibility studies [24]. In the study area, it varies from 0 • to 77 • , and it is a widely acknowledged as being the predisposing factor for debris flows and shallow landslides.
The curvature ( Figure A5) is formed by the intersection of a random plane with the terrain surface [58], and controls the superficial and subsurface water regime of the slopes and, consequently, the predisposition to landslides.
The profile curvature ( Figure A6) is the curvature in the vertical plane parallel to the slope direction. It is the measure of the rate of change of slope and may be planar, concave or convex. Negative values indicate concave portions of the terrain, while positive values indicate convex portions. Terrain with concave curvature are more prone to landslides [59]. A study conducted by Carvalho et al. [60] in a basin of Nova Friburgo municipality showed that the landslides occurred mainly in concave slopes Planar curvature ( Figure A7) is described as the curvature of a contour line formed by intersection of a horizontal plane with the surface [12], and is used to characterize areas of convergence and divergence of surface and subsurface water flow. It is related to the rate of variation of the slope along the contour lines and measures the propensity of the surface and sub-surface flow water to converge or diverge. Positive values correspond to divergent terrain, negative to convergent terrain and null to planar terrain. Convergent or V-shaped profiles are more conducive to trigger mass movements [61].
The flow accumulation ( Figure A8) is derived from a flow direction map, and performs a cumulative count of the number of pixels that naturally drain into outlets and measures the area of a watershed that contributes runoff to each pixel. Flow accumulation highlights the sites that tend to concentrate water after precipitation, which are likely to have a high incidence of landslide triggers [19,62].
Flow accumulation and slope gradient were combined to derive two widely established hydrological indexes, namely the topographic wetness index (TWI)-Equation (1), developed from Beven and Kirkby [63], and the stream power index (SPI)-Equation (2), created by Bagnold [64], which are related to surface and sub-surface water and soil moisture, and thus involved in landslide triggering ( Figures A9 and A10, respectively).
where A = upstream catchment area (m 2 ) or flow accumulation, b is the width of a cell through which water flows (m) and β is the slope angle (radian). High values of TWI are associated with flat sites and valley bottoms, with a large contributing area. These areas are usually related to the deposit zones, with material coming from past landslides. The places where landslides started normally have low TWI values, since they are concentrated in areas with steep slope, also associated with the presence of a thin layer of soil [65]. The SPI is derived from slope and hydrological contributing area. The SPI is higher near the streamlines, and the maximum value occurs where the contribution area is larger, and the morphology is flatter. Then, it is expected that higher SPI values indicate places where runoff is high, which makes it difficult to establish well developed and thick soil layers [65].

Results
The landslide susceptibility indices calculated in each basin were imported in ArcGIS software and used to draw the reclassified susceptibility maps ( Figure 6). The reclassification is based on equal intervals (0-25 low; 25-50 moderate; 50-75 high; 75-100 very high), because this classification scheme is independent from the statistical distribution of the mapped values and, thus, the same class-break values can be used for every map. As a consequence, the differences encountered depend only on the model results, and are not influenced by the classification approach, allowing for a more direct and straightforward comparison among the results of the four models. All the model configurations identify some large susceptibility hotspots around Nova Friburgo and Petropolis, while in many basins scattered away from the most intensely landslide affected areas, the spatial prediction of susceptibility is modelled with different outcomes. Figure 7 can be used to get an idea of the discrepancies between the resulting maps, as it represents the difference between the maximum and the minimum susceptibility values calculated in each basin by the four models. As can be seen, relevant differences can be encountered in some basins. This outcome highlights the sensitivity of the LSM results to the approach used to transform numerical parameters into basin-level variables. The subsequent step of the research consists of assessing which one of the resulting maps is the most reliable.
According to a consolidated literature [24,66], the area under the receiver operating characteristic (ROC) curve (AUC) was used to validate and evaluate the accuracy of each model's results. The ROC curve is a cutoff-independent evaluation model, and the AUC value indicates the performance of the corresponding model configuration, as shown in Figure 8.
The validation statistics show that the most accurate spatial prediction is obtained using the mean value of all numerical variables to characterize each basin (AUC = 74.3%). However, the fourth configuration, which uses all the parameters (maximum values, minimum values and standard deviation), shows a comparable accuracy (AUC = 73.7%). Both outcomes can be considered indexes of a "good" prediction accuracy, according to the scale proposed by Arabameri et al. [67].
As a further test of the modeling results, Table 2 provides a summary of the subdivision of the basins into the four susceptibility classes. These statistics could be used to interpret the meaning of the AUC values derived by the validation of the LSMs. In all four applications, higher susceptibility classes are a lower number, and most of the landslides concentrate there. This is a reasonable outcome for LSM, but some differences can be identified among the four approaches. The LSM derived by the MEAN approach provides the highest AUC, because it is the one with the most balanced classification: it has one of the highest numbers of high and very high susceptible basins, similar to MAXIMUM, but the percentage of correctly classified basin is higher. The use of the MAXIMUM approach leads to an overestimation of hazard (which is reflected by a lower AUC value), as the derived LSM is the one with the highest number of very high susceptible basins and a relevant number of highly susceptible ones, but the number of basins that have actually been affected by landslide are the lowest percentage of all the four approaches. On the contrary, STANDARD DEVIATION underestimate hazard, as the number of basins mapped as highly susceptible is lower than all other approaches.    Another issue worth discussing is the importance of the explanatory variables used in the modeling. The relative importance of each variable can be quantitatively assessed with the "out of bag error" [19,68,69], an index measuring the relative error that would be committed by the model in case an independent variable is not used. Figure 9 shows the importance of the parameters used in the best performing configuration (based on mean values (Figure 9a)) and in the fourth configurations (all parameters (Figure 9b)). The latter can be used to get an overview on how the relative importance of a parameter can increase or decrease if it is aggregated at the basin level considering the mean value, the maximum value or the standard deviation. Indeed, in the first ranks, the parameters of various typology are present: some based on mean values (e.g., mean curvature), some based on maximum values (e.g., maximum slope), some based on standard deviation (e.g., standard deviation of curvature) and some categorical parameters (e.g., soil type).  Table 1 for details on the variables used.

Discussion
A closer investigation of the results obtained can help getting some insights on the use of numerical variables with spatial units larger than the cell size of the raster DEM used to derive them. Indeed, the results showed that the method used to aggregate the values of the parameters over the spatial units of analysis had a noticeable impact on the results, both in terms of overall predictive capability (AUC values) and in terms of spatial distribution of the susceptibility values (maps in Figure 6).
Using the mean value to characterize the whole basin is the strategy that produced the most accurate predictions in terms of AUC (74.3%, as shown in Figure 8). This means that when the environmental information needs to be upscaled from the pixels to relatively larger spatial units like basins, using the averaged values of the hydrologic and morphologic parameters derived from raster DEMs seems to better reflect the average spatial probability of landslide occurrence. In contrast, the model configuration using maximum values to characterize the whole unit basins is the one that returns the lowest AUC values (70.4%).
Indeed, it should be stressed that the approach of using hydrographic basins as a spatial unit is based on the characterization of each basin with a relative and averaged spatial probability of being interested by landslides. Consequently, there is not a full correspondence with the physical mechanism of shallow landslide triggering, for which extreme values are more significant as predisposing factors (for instance, other factors being equal, a shallow landslide is expected more likely when extreme values of morphological parameters are present). In addition, the use of the maximum value alone may not be entirely representative of the geomorphological system of the hillslopes at catchment scale, and could provide biased predictions, in which the hazard is overestimated, hence, the largest number of basins mapped in the highest susceptibility classes (Table 2) and the lowest AUC value obtained after the validation.
Using only the standard deviation derived variables, the overall results are intermediate (AUC = 71.8%). This outcome was quite surprising, as the standard deviation of geomorphic parameters is rarely used in LSM studies requiring spatial aggregations; however, some authors used standard deviation as a proxy for the variability of the conditions encountered in a large spatial unit [19,70,71], and similar indexes may be reasonable indicators of predisposing conditions. Beyond the generalizations inferred by the comparison of AUC values, it should be stressed that the above considerations depend largely on the morphometric parameters considered. For some indexes and attributes with a hydrologic meaning (SPI, flow accumulation and plan curvature) it can be seen from Figure 9b that the parameter derived from the maximum value has a higher importance in the modelling than the homologue parameters derived from the mean or the standard deviation. For elevation, curvature and profile curvature, mean values seems the best method of aggregation (Figure 9a), probably because these parameters are more suited to account for the general characteristics of each basin in terms of morphology or climatic setting; for instance, elevation is clearly related to the mean rainfall amounts falling on the test area and, thus, its average value appears to be a good indicator in the susceptibility analysis. According to the ranking of the importance of the single parameters, slope and TWI are the only factors for which the standard deviation-derived parameter has the highest statistical importance. This outcome can be interpreted as observing that the standard deviation of TWI and slope could be symptomatic of abrupt changes in the hillslope morphology and sub-surface flow. It is well known that the mechanism of triggering in some circumstances can be favored neither by exceptionally high nor by mean values, but by abrupt changes in the hillslope system, such as a marked contrast in geotechnical properties, or sudden changes in surface or subsurface morphologies that reduce hydraulic transmissivity, thus fostering return flows or overland flows that may evolve in soil slips or mobilization of debris. In the case of many other factors, however, the standard deviation is not the most important derivate parameter, but it can still be conveniently used (thus, the second highest performance of the "all-parameters" configuration in terms of AUC) as an index of how much the morphology is varied in each basin, thus giving an approximate idea of the terrain roughness and energy of relief. Indeed, the standard deviation of several morphometric attributes has been already used in landslide susceptibility assessment with this connotation [19,72].
It is not surprising, then, that a good spatial prediction of landslide susceptibility is obtained with the configuration, encompassing all the factors derived from the maximum, mean and standard deviation (AUC 73.7%); this strategy exploits each approach used to aggregate pixel values at the basin level, and can also be considered a further confirmation of the flexibility and predictive power of random forest, which can handle very complex spatial information, even when provided by a large number of input variables affected by mutual relationships. Among the 28 used parameters, some resulted to have little importance; however, their contribution was still evaluated positively by the internal tests of the RF algorithm (Figure 9b). To further corroborate this assumption, some model configurations encompassing a reduced set of parameters have been tested and validated: the validation statistics showed that a more parsimonious parameterization does not improve the model performance: the AUC value is reduced to 70.3% if the least important explanatory variable (Flow_Acc_STD) is excluded from the analysis, and to 68.8% if the second least important (Flow_Acc_MEAN) is also excluded. Similarly, a configuration in which the redundancy and the parameterization is reduced to the minimum was tested: for every parameter only the variable that showed the highest OOBE value was selected. This configuration therefore uses each parameter only once, and combines the variables aggregated using the mean value (elevation, curvature and profile curvature), some using the maximum (planar curvature, SPI and flow accumulation) and some using the standard deviation (slope and TWI) and gets lower AUC scores (0.700) than the four main configurations tested. However, the AUC value reported by the "ALL" configuration is lower than the "MEAN" configuration, suggesting that the random forest algorithm is not completely capable of resolving the problems of collinearity between such a high number of variables.
Another important input for discussion is the quality of the results obtained. The validation of "ALL" and "MEAN" configurations returned the AUC scores 74.3% and 73.7%, respectively. These values fall in the 70-80% range, which is considered "good", according to the classification proposed by some authors [67]; however, it is lower than validation values commonly reported in the international literature. This indicates the limits of the proposed approach, in which the predisposing factors at the triggering point of each landslide cannot be ascertained and included in the analysis, because only the general characteristics of each basin is taken into account. This approximation has been discussed before, here, it is worth stressing that some predisposing factors dramatically lose their significance if considered according to this approach: for instance, the flow accumulation has a very low predictive importance for every aggregation method (as reported in Figure 9b, Flow_Acc_MEAN and Flow_Acc_STD are the least important variables, and even Flow_Acc_MAX has a very low relative importance).
In addition, another reason for the relatively low AUC scores is that the use of basins as spatial computation unit reduces the training and validation sample. In the study area the landslide inventory contains a high number of landslides (namely, 3483). In a pixel-based approach, that would allow to identify at least 3483 points to sample predisposing conditions (of course, only if all the landslide locations were known with good accuracy). In the proposed approach, only 1253 basins were defined, and only 312 of them included landslides; consequently, the training and validation samples are smaller. It has not been excluded that in the future, with a larger and updated dataset, the results could be improved. However, it should be noted that a pixel-based approach could not be pursued in this test site, due to limitations in the landslide dataset: even small uncertainties in the locations and geometry of the landslides prevent all of them from being related exactly to a pixel, and therefore a pixel-based analysis would be impossible.
Another main drawback of the proposed approach is the coarse spatial resolution represented by the hydrographic basins. However, this choice was made for the practical reasons of integration in the framework of existing and planned hazard management procedures, since it is a good practice to consider also the needs of decision makers in the planning of a research work [73]. Raster cells would allow for a finer spatial resolution, but cannot be used as a spatial unit for warnings and operative response; thus, a spatial reaggregation in larger spatial units is needed. Basins would be easier to understand for populations, and untrained personnel and could be conveniently exploited in decision making and in operational risk management procedures.
Indeed, the strict collaboration between CEMADEN and local civil defense authorities was useful to direct the research towards practical and operative solutions. As a result, the final LSM will be used as an operational tool for hazard management. Regarding risk prevention activities, the identification of basins with higher landslide susceptibility will be useful in terms of setting the priorities of intervention for risk mitigation strategies. Use of the LSM is also planned during the prediction and hazard response activities to be undertaken in the situation room: joint use with weather monitoring and forecasting products will allow for a better constraint of the spatial extent of the incoming hazardous events. However, regarding this operational use, it should be stressed that, like any other LSM, the resulting map is a static instrument that provides a spatial assessment of the relative spatial probability of landslide occurrence, based on the analysis of predisposing factors. Therefore, triggering factors (e.g., rainfall intensity) are not directly accounted for, and no temporal information is provided. To overcome this drawback common to all LSMs, in the future, progress of the research and integration of the susceptibility map with a warning system based on rainfall thresholds [40] is planned, using the recently proposed combination techniques [73][74][75].
Moreover, the option of basin-scale aggregation is the only option that makes it possible to handle a multi-source landslide inventory with inhomogeneous characteristics (e.g., different landslide geometries and mapping approaches). If all landslides are associated with the hydrographic basin spatial unit, it makes no differences if landslides are mapped as points or polygons, and if the element mapped is the scar, the crown or the whole accumulation body, as long as, of course, each landslide begins, evolves and ends within the same catchment. In addition, small uncertainties in the spatial location of the mapped elements also become irrelevant, or at least would have a much lower influence than in a pixel-based approach.

Conclusions
In this study, the random forest technique was applied to map landslide susceptibility in a mountain region north of Rio de Janeiro, Brazil, where hundreds of landslides occurred in 2011, and impacted the environment and the population, causing hundred deaths and economics losses.
The study area was divided into drainage basins to be used as meso-scale spatial unit for the analysis and the resulting map. This approach was chosen for two main reasons. First, the choice of this spatial unit is mainly linked to the operational hazard management procedures in which the resulting map is planned to be used: basins are coarse spatial units, but they have the advantage of referring to a portion of the territory that is easy to understand, and in which a unitary counteraction can be implemented by Civil Defense organizations in case of alarm. Second, the hydrographic basin approach makes it possible to overcome the problems of spatial accuracy and inhomogeneity encountered in the multi-source landslide inventory used; such limitations would prevent a pixel-based analysis, because the uncertainties in locations and geometry of the landslides would make it impossible to exactly relate each landslide to a pixel.
Twenty-eight variables, derived from four categorical and eight numerical parameters (described in Table 1), have been used for the model construction. To aggregate their values from the pixel to the basin level, different approaches were pursued. For categorical parameters, the most representative (majority) class was assigned to each basin, and for the numeric parameters, three different aggregation schemes were implemented, based on the maximum value, mean value and the standard deviation. Four susceptibility maps were generated varying the approach used to characterize the numerical morphometric parameters: maximum values, mean values, standard deviation and all values together. The best results were obtained using mean values (AUC 74.3%) and all values together (AUC 73.7%), while the standard deviation and maximum values produced lower validation scores (AUC was 71.8% and 70.4%, respectively).
The LSM obtained with this work will be used by CEMADEN and local civil defense authorities as an operational tool for hazard management. The identification of basins with higher landslide susceptibility will be useful to set priorities of intervention for risk mitigation strategies. Moreover, joint use in the situation room with weather monitoring and forecasting products will allow for a better constraint of the spatial extent of the incoming hazardous events.
In addition to producing a meso-scale susceptibility map for the study area, the present work made it possible to draw some conclusions concerning a common source of uncertainty that affects LSM based on spatial units larger than the pixel size of the input parameters: that is, how to generalize the fine-level pixel information at the level of a larger spatial unit like basins:

•
In the analyzed case of study, with peculiar physiographic characteristics and a limited number of training data, the best generalization approach was the use of the mean value of each numerical parameter to derive the explanatory variable to be used in the susceptibility analysis. • Nevertheless, investigating the relative importance of each variable, it was discovered that each morphometric factor could be more conveniently exploited using a different generalization scheme, related to the different physical meaning of each of them.

•
For instance, in the proposed application, most of the hydrologic factors are better represented by the maximum value, elevation and curvature by the mean value, slope and TWI by standard deviation.

•
Given that, if the susceptibility assessment is based on a model that can handle multicollinearity, large number of parameters and complex configurations, another option worth testing could be using all the parameters derived with the three different approaches.