Development of a Landslide Hazard Index to Land Use by Using Logistic Regression and Geographic Information System: A Case Study in Pingtung County in Southern Taiwan

1. Introduction

Taiwan is an island located in the North Pacific subtropical monsoon region, which is vulnerable to typhoons. More than two-thirds of Taiwan’s total area is mountainous with steep slopes. Recently, climate anomalies have resulted in frequent extreme rainfall events. The number of typhoons that hit Taiwan has increased every year, and large-scale landslides and debris flows have become likely. As stated in the 2012 White Paper on Disaster Management [1], remote mountainous settlements in Chiayi County, Kaohsiung City, Tainan City, Pingtung County in southern Taiwan account for approximately 24% of the total area in Taiwan. Recent large-scale sediment disasters due to typhoon-induced concentrated rainfall occurred mostly in the remote mountainous areas of southern Taiwan. Such disasters pose a major threat to the environmental security of the settlements in remote mountainous areas. Since Typhoon Morakot struck in 2009, the slopes of the remote mountainous areas of southern Taiwan have become highly unstable. Therefore, studies on the land use risk of the aforementioned settlements during rainfall disasters are essential.

Jan and Lee [2] used the product of effective accumulated rainfall (EAR) and rainfall intensity as a rainfall trigger index (I_RT) for determining debris flow; they discovered relevant geological factors (average stream bed slope, catchment area, landslide rate, and lithology) for establishing comprehensive geological indicators. A rainfall prediction model for debris flow was constructed using the rainfall warning and comprehensive physiographic indices. For rainfall warning, the following occurrence possibilities were set: 10%, 50%, and 90%. Chen et al. [3] demonstrated that extreme rainfall events causes frequent landslide erosion; they reported that a maximum 24-h rainfall amount of >600 mm resulted in the average landslide erosion rate of 64%–79%. Jan et al. [4] mentioned that the effects of typhoon rainfall depend not only on the rainfall amount but also on its intensity. Tseng et al. [5] reported that the number and area of rainfall-induced landslides are positively correlated with the degree of land disturbance. Chen et al. [6] investigated areas with high and low susceptibility to landslides; they identified cumulative rainfall to be the primary factor for landslide occurrence in areas with high susceptibility; by contrast, in areas with low susceptibility, the rainfall intensity was found to be the key factor. Because of climate extremes, areas with high susceptibility exhibited a higher magnitude and frequency of landslides than did those with low susceptibility.

Lee et al. [7] explored the mechanisms underlying landslide disasters due to the intensity of Typhoon Soudelor–induced rainfall on the mountainous areas of northern Taiwan; they revealed that most sediment disasters occurred 1 or 2 h after peak rainfall. Moreover, the maximum rainfall occurring 3, 6, and 12 h after this event was higher than that in the 200-year return period. To estimate a rainfall alert value for landslide occurrence, Caracciolo et al. [8] analyzed historical data on the rainfall that led to a landslide as well as data on rainfall occurring 5, 15, and 30 days before a landslide in southern Sicily, Italy. Chen and Wu [9] focused on a state-owned forest land in Taiwan. They analyzed various internal (e.g., elevation, slope, terrain roughness, distance from the fault, and distance from the river) and external (total accumulated rainfall after event) factors through logistic regression (LR) to build an internal potential level without external factors; subsequently, they evaluated the changes in the potential value and rainfall to estimate a reference value for the prediction of landslides induced by rainfall. Shahabi et al. [10] combined remote sensing and geographic information system (GIS) to statistically delineate areas susceptible to landslides; for this, they considered the following factors: slope, aspect, elevation, lithology, normalized difference vegetation index, vegetation, rainfall, distance from the fault, distance from the river, and distance from the road. Tseng et al. [5] focused on the periods before and after a typhoon that recently caused a road slope landslide in the study area; through image interpretation, they identified surface changes before and after the landslide disaster and developed a model for assessing landslide susceptibility. Then, the GIS platform was used to construct landslide susceptibility maps. Among the statistical methods used for evaluating landslide susceptibility, LR analysis has been proven to be one of the most reliable approaches [11,12,13,14,15,16,17,18].

In the present study, we used the random tree classifier (RTC) to interpret and classify the images obtained for the periods before and after rainfall. The findings can provide information on surface land use and disasters, which can be used to identify the areas susceptible to sediment disasters in Sandimen and Wutai townships, Pingtung County. Considering the regional environmental characteristics and land disturbance, we further analyzed the potential of slope environmental intensity, factors associated with slope disasters, and rainfall characteristics of the study area and investigated the correlation between rainfall trigger and slope environmental strength potential to develop I_RL. Our findings ensure the safe use of land in mountainous settlements.

2. Methods

2.1. Rainfall Projection

2.1.1. EAR

We referred to the rainfall analysis methods described by Seo and Funasaki [19] and Tseng et al. [5]. Concentrated rainfall is considered to be continuous if no rainfall occurred 24 h before and after the rainfall event. The rainfall field with continuous rainfall inducing a landslide is regarded as the main rainfall field. The beginning of rainfall is defined as the time point when the first rainfall reached ≥4 mm in the main rainfall field. Cumulative rainfall inducing a landslide is calculated. Cumulative rainfall is divided into previous indirect rainfall (P_b) and previous direct rainfall (P_r). Previous indirect rainfall refers to the amount of rainfall in the main rainfall field within 7 days [19]. It can be calculated using Equation (1).

$${\displaystyle {\sum }_{n=1}^7{k}^n{P}_n=P_b}$$

(1)

P_n is the amount of rainfall (mm) in n days before the main rainfall field; k is the decreasing coefficient. In the present study, k is 0.9 [20]. EAR can be calculated as follows:

$$EAR=P_r+P_b$$

(2)

P_r is the amount of cumulative rainfall from the first rainfall in the main rainfall field from the time of the landslide disaster (a time point when the first rainfall reaches ≥4 mm) to the occurrence of the landslide.

2.1.2. Rolling Rainfall Intensity

Rain-induced landslides may be triggered by continuous rainfall for several hours. Therefore, rolling rainfall intensity can be expressed using Equation (3):

$$I_{mR}={\displaystyle {\sum }_{i=t-m+1}^m{I}_i}=I_{t-m+1}+I_{t-m+2}+\cdots +I_m$$

(3)

I is the rainfall intensity, m is the unit time of rain rolling, and m = 3 h [5]. I_mR is the I_R in m h, and I_t is the rainfall intensity in t h.

2.2. RTC

We used the RTC in the ArcGIS supervised image classification module to classify images. The RTC does not lead to overfitting and can process segmented images and auxiliary grid data sets [21]. The classifier can be used to construct several decision trees, and it selects a random subset of variables for each tree and uses the most frequent tree output for the overall classification. Therefore, random trees correct for the tendency of decision trees to overfit training data. Random trees are a collection of individual decision trees; each tree is generated from different subsets of sample and training data. Decisions are made in the order of the importance of each pixel being classified. An image is drawn for a certain pixel appears as a tree branch. After the entire data set is processed, each branch forms a whole tree; this is the concept of random trees. In the operation of random trees, each tree has decision-making power; this process reduces overfitting. In the random tree method, analogically, many trees continue to grow; the changes in the trees of a forest are projected into randomly selected subspaces through training data, and each node’s decision is optimized through a random process [21].

2.3. LR

In regression analysis, dependent variables serve as categorical variables, whereas independent variables serve as continuous or dummy variables. Dummy variables represent categorical data; the corresponding numerical values are the basis for classification and have no comparative significance. The primary feature of LR is that dependent variables are categorical variables, whereas independent variables can be continuous or categorical variables [22].

LR is used for analyzing dichotomous dependent variables. One or more independent variables may be included in the model. Response variables can be categorical or continuous variables. Previously, landslide susceptibility was analyzed mostly considering the distribution patterns of binary variables [23]. Through a statistical induction method, a set of regression patterns may be identified to differentiate between landslide and nonlandslide. The value range of the logistic distribution function is 0–1, and the distribution of the value range follows an S-shaped curve [24].

As shown in Equation (4), linear regression assumes a linear function and includes dependent and independent variables (and random residual values) [22].

$$Y={\beta }_0+{\beta }_1\times x_1+{\beta }_2\times x_2+{\beta }_3\times x_3+e$$

(4)

Y is a dependent variable; x₁, x₂, and x₃ are independent variables; and e represents a residual value. β₀+β₁x₁+β₂x₂+β₃x₃ is a linear function between expectation E(Y) and the three independent variables.

LR is used to construct practical and reasonable allocation models for predicting the correlations between a dependent variable (y) and a set of independent variables (x). The relationship usually refers to the relationship between a set of independent variables (x_s), which is used to predict the probability of the dependent variable being equal to 1, such as the probability of landslide occurrence. The ratio of the probability of event occurrence to that of event nonoccurrence is called the event odds. A linear function can be obtained using its natural logarithm, which is the logit model, as shown in Equation (5) [22].

$$\ln \left(\frac{P}{1-P}\right)=\alpha +{\displaystyle \sum _{i=1}^k{\beta }_i{x}_i}$$

(5)

P is an independent variable; x₁, x₂,….x_k are the probabilities of event occurrence. In this study, y indicates the probability of landslides, and the aforementioned x parameters represent various independent variables. A logistic curve can be constructed using the logit function; its mathematical formula is shown in Equation (6).

$$P=\frac{1}{1+e^{-(\alpha +{\displaystyle \sum _{i=1}^k{\beta }_i{x}_i})}}$$

(6)

In general, the probability threshold is set at .5. A predicted probability (p) of ≥.5 indicates the likelihood of a landslide event; by contrast, if p < 0.5, a landslide event is not expected to occur.

3. Study Area

This study was conducted in Sandimen Township and certain mountainous areas of Wutai Township in Pingtung County in southern Taiwan (Figure 1). Pingtung County is located in the south of the Tropic of Cancer. Except for elevated mountainous areas, the county has a tropical monsoon climate. The temperature does not vary substantially throughout the year. The annual average temperature is approximately 25.5°C. It is humid and rainy in summer. However, because of the barrier between the Dawu Mountain Range and the Central Mountain Range, the cold northeast monsoon is blocked from entering the county; this blockage combined with a low latitude and sufficient sunshine reduces the intensity of winter in Pingtung, where the average high and low temperatures are 24–27°C and 16–19°C, respectively [25,26]. In Sandimen and Wutai townships, most of the strata are Chaozhou and Bilushan formations; low contents of alluvial and grounding deposits are observed. Lithologically, the study area is mainly interbedded with hard shale, slate, and sandy shale. In addition, low contents of gravel, sand, clay, and slate interbedded with metamorphic sandstone and metamorphic igneous rock lens, mud, sand, and gravel, lenticular sandstone bodies are noted [27]. Sandimen Township is located at an altitude of 100–2,159 m and at the intersection of mountains and plains and features mostly hilly terrain. Wutai Township is located in the northeast area of Pingtung County (in the Central Mountain Range). This township has steep terrain and is located at an average altitude of >1,000 m. The largest tributary flowing through the study area is Ailiao River, which is the largest tributary of Gaoping River. Its farthest source stream is Eluowu River, which flows south to Babanaban River and then turns southwest to Laibuan River (where it is called Ailiao North Stream); subsequently, it turns south to join Hayou Creek and Qiaoguo Laci Creek (among others) and finally turns south to converge into Ailiao South Stream in Dalai Village; it is called Ailiao Creek in Sandimen Township. The main rivers in Wutai Township are Erchong and Ailiao Beixi, whereas the main rivers in Wutai Township are Ailiao Beixi and Qiaoguolaci [25].

In Sandimen Township, Sandi and Dalai villages have the highest (n = 1,683) and lowest (n = 402) populations. In Wutai Township, Wutai and Jilu villages have the highest (n = 1,325) and lowest (n = 210) populations [28]. The residents of Sandimen Township are mainly Paiwan individuals (Taiwan aborigines); some Rukai individuals also reside in this township, such as in Qingye Village, which has recently been included in the Maolin National Scenic Area. The main economy in Sandimen Township is agriculture, and the main agricultural products are sweet potato, millet, taro, mango, pineapple, red quinoa, and coffee [29]. The residents of Wutai Township are mainly Rukai individuals [30]. The main agricultural products of the township are love jade, red quinoa, millet, coffee, taro, and sweet potato [31].

4. Image Interpretation and Classification

4.1. Preprocessing of Satellite Images

FORMOSA 2 (FS-2) or SPOT-5 satellite images of the areas affected by typhoons or rainstorms were acquired. Data on the 2009 Typhoon Morakot, 0517 rainfall in 2013, 2013 Typhoon Soulik, and 2013 Typhoon Kongrey were obtained to explore the land development types and landslide lands in the study area. We used six satellite images with image resolutions of 8 m × 8 m and 10 m × 10 m (Table 1).

The FS-2 and SPOT-5 satellite images are multispectral and contain four spectral bands [32]. FS-2 spectral bands are red, green, blue, and near-infrared (NIR), whereas SPOT-5 spectral bands are red, green, NIR, and short-wave infrared. The telemetry image processing software ERDAS IMAGINE [33] was used to fuse and locate the images.

4.2. Selection of Satellite Imagery Classifications

In this study, FS-2 or SPOT-5 satellite images of four typhoons or rainstorms were selected as the base map. Using GIS ArcGIS, each classification factor was manually digitized and circled. We compared the Tiff files of the satellite images (resolution, 2 m × 2 m) of the study area to increase the accuracy of sampling area delineation; furthermore, we selected the following eight classifications suitable for the scope of this study: river, road, building, farmland, forest, meadow, streamway, and bare land (Figure 2).

4.3. Interpretation of Images and Assessment of Accuracy

Considering the aforementioned eight classifications, the RTC module was used with the results of texture analysis to interpret the images. The research scope was interpreted using a total of six satellite images obtained before and after the 2009 Typhoon Morakot, 0517 rainfall in 2013, 2013 Typhoon Soulik, and 2013 Typhoon Kongrey. Figure 3 depicts the results of image interpretation. To confirm the accuracy of the interpretation results, we randomly selected a total of 25 points (interpretation grids; checkpoints) for each classification and compared them with high-resolution aerial photos and on-site survey data. The most common accuracy evaluation method is the error matrix [34], which is used to calculate the coefficient of agreement (Kappa index) and overall accuracy. The value of the Kappa index proposed by Cohen [35] ranges from 0 to 1. Kappa values of <0.4, 0.4–0.8, and >0.8 indicate low, medium, and high accuracies, respectively. Table 2 summarizes OAs before and after the rainfall events in the study area. In this study, the average Kappa value of each satellite image was approximately 0.71, and the average overall accuracy was approximately 74% (medium to high accuracy).

4.4. Identification of Landslides Through Image Interpretation

After image classification, to identify rainfall-induced landslides in the study area, we subtracted bare land grids before and after rainfall events by using the image subtraction method. The river, streamway, and bare land classifications with a slope percentage of <5% were deducted. High-resolution aerial photos of the study area were compared, and manual inspection was performed to identify landslides; then, the locations of the rainfall-induced landslides in the study area were obtained (Figure 4). The 2009 Typhoon Morakot resulted in the highest rainfall and the largest landslide area (1,313.12 ha). The areas of landslides induced by the other events were as follows: 0517 rainfall in 2013, 813.92 ha; Typhoon Kongrey, 789.44 ha; and Typhoon Soulik, 635.04 ha).

5. Development of a Model for Assessing Rainfall-Induced Landslide Susceptibility

5.1. Selection of Factors Associated With Landslide

Landslides can be attributed to factors such as artificial slope land use, environment, and rainfall triggers. Human activities may negatively affect water and soil conservation on slope lands, which compromises the safety of slope land. Therefore, referring to studies of Chen et al. [36] and Tseng et al. [5], we divided the aforementioned factors into the following three categories: environmental, slope disturbance, and rainfall trigger factors. The environmental factors assessed in this study included elevation, slope, aspect, slope roughness, terrain roughness, distance from the river, and geology. The slope disturbance factors included road density, building density, farmland planting rate, forest density, grassland density, and bare density. The rainfall trigger factor was the product of EAR and the maximum 3-h rolling rainfall intensity (I_3R,max).

The digital elevation model (DEM) was used to analyze the environmental factors using ArcGIS Spatial Analyst. A basic grid of 40 m × 40 m was constructed for the factors related to landslide occurrence, followed by the development of a GIS database of the study area. The hazard factors are described below.

5.1.1. Environmental Factors

5.1.1.1. Elevation

The elevation of a place refers to its height from the sea level. In Taiwan, the elevation datum comprises the mean sea level of Keelung Port, and the elevation of the terrain is called the elevation difference. When other factors are constant, higher degrees of elevation are associated with a higher sliding force. We used a 20 m × 20 m DEM for analysis. The elevation of the study area is 78–2,437 m, which was divided into a total of seven grades and coded. Table 3 presents the elevation classification codes, ranging from an elevation of <350 m coded as 1 to an elevation of >2,100 m coded as 7.

5.1.1.2. Slope

Here, slope indicates the inclination of a slope. Higher degrees of inclination indicate steeper slopes and vice versa. High inclination as well as poor soil and water conservation may facilitate the development of landslides. In this study, the slope value of each grid was obtained through DEM and ArcGIS slope analysis. According to the Classification Standards for Land Utilization Limits of Hillside Lands outlined in the Regulations of Soil and Water Conservation [37], we graded and coded slopes (Table 4). A first-grade slope was coded as 1, whereas a seventh-grade slope was coded as 7.

5.1.1.3. Slope Roughness

Slope roughness reflects the changes in a slope (standard deviation). Higher degrees of slope roughness indicate higher extents of changes, which may facilitate the occurrence of landslides. Cluster analysis was performed to classify the slope roughness of the study area. Slope roughness of <18.15 was coded as 1, whereas that of >63.34 was coded as 7 (Table 5).

5.1.1.4. Aspect

Each slope has a different aspect. Aspect refers to the inclination direction of a slope. It affects wind flow and rainfall distribution. In this study, using DEM and ArcGIS aspect analysis, the aspect data of each grid were obtained. The inclination angle may be as follows: north, northeast, east, southeast, south, southwest, west, northwest, and flat ground (clockwise). In this study, a flat ground was coded as 1, whereas the southwestern part of a windward side was coded as 6 (Table 6).

5.1.1.5. Surface Roughness

Surface roughness reflects the changes in a surface (standard deviation). Higher degrees of surface roughness indicate higher extents of changes, which may facilitate landslide occurrence. In the present study, cluster analysis was performed to classify the surface roughness of the study area. Surface roughness of <417.62 was coded as 1, whereas that of >2,298.55 was coded as 7 (Table 7).

5.1.1.6. Distance From the River

The distance of an area from a water body may affect the occurrence of landslides. A shorter distance is associated with higher levels of groundwater, which softens soil, thereby facilitating landslide occurrence. We used a river map to calculate the distance of each grid from a river; for this, the ArcGIS buffer analysis function was used. The analysis results were coded as shown in Table 8.

5.1.1.7. Geology

We used the geological map of the Central Geological Survey of the Ministry of Economic Affairs (2021) to obtain information on the geological conditions corresponding to each grid in the study area. On the basis of geological age and lithological data, it was divided into four strengths: very weak, medium, medium strong, and strong [38]. The highest strength was coded as 1, whereas the lowest strength was coded as 4 [36].

5.1.2. Slope Disturbance Factors

Using RTC with texture analysis for image interpretation and classification, we assessed the slope disturbance of the research area. The following six factors were included in further analysis: road density, building density, farmland planting rate, forest density, grassland density, and bare density. To quantify the slope disturbance in each basic grid, slope disturbance was defined as the area percentage ratio of each disturbance factor in each basic grid. The degree of land development of hillsides and its index were evaluated in reference to earlier studies [5,36]. Because we regarded environmental and slope disturbance factors as independent factors influencing landslides, a slope disturbance index was developed and is shown in Equation (7).

$$I_{DC}={\displaystyle \sum G_{DC}}\times R$$

(7)

DC refers to the disturbance condition. G_DC is a DC grade, which is the score corresponding to slope land use disturbance in each grid. R is the ratio of the area occupied by each slope and the utilization factor corresponding to the grid. Table 9 presents the grades corresponding to the DCs of the aforementioned six slope disturbance factors. For grading, we referred to a study conducted by Chen et al. [36].

5.1.3. Rainfall Trigger Factors

We obtained relevant data from the following 12 rainfall measurement stations within and near the study area: Weiliaoshan, Majia, Hongyeshan, Shangdewun, Dajin, Jinfeng, Ligang, Meinong, and Gusia [26] as well as Zhiben-5, Ali, and Sandimen [39]. Through inverse distance weighting (ArcGIS), the EAR and I_3R,max values were calculated for each station (Table 10 and Table 11) for four rainfall events occurring in 2009 and 2013.

5.2. Analysis of Slope Environmental Strength Potential

We performed correlation analyses of environmental and slope disturbance factors. After highly correlated factors were eliminated, LR was performed to analyze slope environmental strength; thus, a slope environmental strength regression model was constructed. Then, we established the index of slope environmental strength potential (I_SESP). A basic grid of 40 m × 40 m was constructed using ArcGIS; subsequently, the slope environmental strength potential value in each grid of the study area was calculated.

5.2.1. Correlation Analyses of Slope Environmental Strength Potential Factors

The slope environmental strength potential factors initially selected in this study included environmental factors such as elevation, slope, aspect, slope roughness, surface roughness, distance from the river, and geology and slope disturbance factors. To ensure that the factors are independent and exert no mutual effects, we performed correlation analyses of the selected factors before LR analysis to identify the degree of correlations between the factors. Elevation, slope, slope roughness, surface roughness, and slope disturbance index served as continuous variables, whereas aspect, distance from the river, and geology served as categorical variables. In LR, categorical variables only represent the data distribution of factor codes and not the effects of the code size on dependent variables. Hence, we focused only on continuous variables. Statistical analysis was performed using SPSS [40]. For the four rainfall events, the correlations between the factors were investigated through Pearson correlation analysis. If the absolute value of the correlation coefficient is closer to 0, the degree of correlation is weaker; by contrast, if it is closer to 1, the degree of correlation is stronger. Table 12 presents the correlations between slope environmental strength potential factors. Elevation exhibited high-level correlations with surface roughness. Slope exhibited medium-to-high-level correlation with slope roughness. The remaining factors exhibited low-level correlations. Hence, surface roughness was not included in the LR model.

5.2.2. Results of the Analysis of Slope Environmental Strength Potential

To reduce errors due to the evaluation model and subjective errors due to selection bias (human errors), we randomly sampled the same number of landslide and nonlandslide samples. After sampling, the data were randomly divided into training (70%) and test (30%) data sets. An environmental strength assessment model was constructed through LR [40].

The dependent variable was landslide or nonlandslide, and the independent variables were environmental factors, such as slope, aspect, elevation, distance from the river, slope roughness, geology, and slope disturbance factors. Through LR, a regression formula was developed for analyzing slope environmental strength potential. The formula is shown in Equation (8).

$$L=\alpha +{\beta }_1\times El+{\beta }_2\times Sl+{\beta }_3\times Slr+{\beta }_4\times I_{DC}+{\displaystyle \sum _{i=1}^4{\beta }_{As,i}{D}_{As,i}}+{\displaystyle \sum _{i=1}^6{\beta }_{Dr,i}{D}_{Dr,i}}+{\displaystyle \sum _{i=1}^3{\beta }_{Gs,i}{D}_{Gs,i}}$$

(8)

El denotes the elevation, Sl is the slope, Slr is the slope roughness, I_DC is the slope disturbance index, As is the aspect grade, Dr is the distance from the river, Gs is the geological grade, β is the regression coefficient of continuous variables, and α is a constant. D represents a dummy variable, which is a value coded for the category corresponding to a factor. Substituting each regression coefficient and constant value into Equation (8) and substituting it into Equation (9), we can calculate the probability value of each grid, which represents the value of I_SESP.

$$P=\frac{1}{1+e^{-L}}$$

(9)

LR was performed to evaluate the training classification results of various rainfall events. The accuracies of classifications performed using training and test data after 2009 Typhoon Morakot were 68.3% and 68.1%, respectively. Training data were used to deduce the accuracies of classifications for the entire study area. The estimated accuracy was 65.2% (Table 13). Table 14 presents the errors in the evaluation of various rainfall events. The accuracies of classifications performed using training data after 0517 rainfall, Typhoon Soulik, and Typhoon Kongrey were 83.8%, 86.2%, and 86.4%, respectively; the corresponding values for test data sets were 84.3%, 86.4%, and 87.5%. The average accuracy of the overall classification was approximately 80.4%.

5.3. Establishment of a Landslide Hazard Index

Higher values of slope environmental strength potential are associated with larger amounts of rainfall and thus higher probabilities of landslides. By referring to the study of Jan and Lee [2], we adjusted it. I_RT can be calculated as shown in Equation (10).

$$I_{RT}=EAR\times I_{3R,\max }$$

(10)

We established a landslide hazard index (I_RL) using I_RT and I_SESP. The formula is shown in Equation (11).

$$I_{RL}=I_{RT}\times I_{SESP}$$

(11)

Higher values of I_RL indicate a higher probability of landslides. At a constant I_RL, I_RT is inversely proportional to I_SESP. Thus, if I_SESP is relatively stable, a large amount of rainfall would induce landslides. By contrast, if I_SESP is low, only a small amount of rainfall would induce landslides. We used the following comprehensive indicators of rainfall-induced landslides [2]: I_RL1, I_RL10, I_RL25, I_RL50, and I_RL90. The index development method is described below.

• I_RL1: From the I_RL values corresponding to the grid data of all rainfall-induced landslides in the study area, the value with a cumulative probability of 1% is selected (Weber’s method) and is indicated as I_RL1. Grids with I_RL values less than that of I_RL1 have landslide probabilities of <1%.
• I_RL10: From the I_RL values corresponding to the grid data of all rainfall-induced landslides in the study area, the value with a cumulative probability of 10% is selected (Weber’s method) and is indicated as I_RL10. Grids with I_RL values between the values of I_RL1 and I_RL10 have landslide probabilities of 1%–10%.
• I_RL90: Landslide and nonlandslide grids with values less than that of I_RL10 are excluded. From the I_RL values corresponding to the remaining grids, the value with a cumulative probability of 90% (Weber’s method) is selected and indicated as I_RL90. Grids with I_RL values exceeding that of I_RL90 have landslide probabilities of >90%.
• I_RL10–I_RL90: To determine a landslide probability corresponding to a comprehensive index between I_RL10 and I_RL90, the relationship between the aforementioned index and landslide probability can be expressed as shown in Equation (12).

$$\frac{I_{RLi}-I_{RL10}}{I_{RL90}-I_{RL10}}=\frac{P-0.1}{0.9-0.1}$$

(12)

Equation (12) can be rewritten as Equation (13).

$$I_{RLi}=I_{RL10}+\Delta I_{RL}\left(\frac{P-0.1}{0.8}\right)$$

(13)

$\Delta I{}_{RL}=I_{RL90}-I_{RL10}$. Thus, the landslide probability corresponding to the landslide comprehensive index I_RLi can be obtained.

5.4. Analysis of Landslide Hazard to Land use

5.4.1. Correlation Between Slope Environmental Strength Potential and Rainfall-Induced Landslide

We first explored the distribution of landslide and nonlandslide grids with the same I_RL value for each event and performed stratified random sampling according to the ratio of the numbers of landslide and nonlandslide grids with the same I_RL value. To calculate the I_RL values corresponding to I_RL1, I_RL10, I_RL25, I_RL50, and I_RL90, after merging the data points of each rainfall event, we sorted the landslide grids; then, using Equation (12), the I_RL values corresponding to different cumulative probabilities were calculated. Table 15 presents the I_RL indices and the corresponding values.

Using the I_RL indices presented in Table 15, the correlation between I_RT and slope environmental strength potential was calculated after merging the data points of each rainfall event (Figure 5). The product of I_RT and I_SESP was the closest to the grids with different cumulative probabilities. The trend line corresponding to each I_RL value was drawn using the power method. As shown in Figure 5, the I_RT of 2009 Typhoon Morakot was the highest among the four rainfall events; thus, most of its landslide grids were distributed above I_RL50. For 0517 rainfall in 2013, the grids were mostly below I_RL90. The I_RT of 2013 Typhoon Soulik was the lowest among the four rainfall events; most of its landslide grids were distributed below I_RL25. The I_RT of 2013 Typhoon Kongrey was higher than those of Typhoon Soulik and 0517 rainfall events; for Typhoon Kongrey, the grids were distributed between I_RL1 and I_RL90.

5.4.2. Estimation of Landslide Hazards to Land use

After the classification of I_RL indices with different cumulative probabilities, the numbers of landslide and nonlandslide grids with different I_SESP grades and the landslide ratio (landslide /nonlandslide) were estimated. Considering that higher values of I_RT are associated with higher values of I_SESP (and higher degrees of hazard), we investigated the most reasonable average distribution of the landslide ratio in an inductive manner. Then, the degree of hazard was calculated for each distribution situation and substituted into the grids of each event. Thus, the hazard degree of rainfall-induced landslide in each grid of each event was obtained. Table 16 presents the landslide ratios (before and after induction) of various I_SESP grades of different I_RL indices. The same color in the table represents the same category of landslide ratio intervals (n = 8). Figure 5 illustrates the results presented in Table 16. The landslide ratio intervals were normalized using the average landslide ratio of the same interval, as shown in Equation (14).

$$Z_{norm}=\left(\frac{X-X_{\min }}{X_{\max }-X_{\min }}\right)$$

(14)

Z_norm represents the value after normalization, X represents the value to be normalized, X_min represents the minimum value in the data, and X_max represents the maximum value in the data. In this study, the normalized value was regarded as the degree of rainfall-induced landslide hazard. The estimated degree of hazard was substituted into the grid of each rainfall event to obtain the degree of hazard corresponding to each rainfall event grid in the study area. The normalized value is a value between 0 and 1. Landslide probabilities of (I_RL) of <1% still represent the possibility of landslide. Therefore, the minimum value of I_RL1 indicates its risk probability. Table 17 presents the hazard value of each landslide ratio interval. The color reference in Table 17 is the same as that in Table 16. The hazard value of each grid of each rainfall event was substituted into each grid, and hazard maps were constructed. Higher degrees of hazard indicate higher probabilities of landslide and vice versa.

The hazard value of each interval was substituted into each grid of each rainfall event, and the hazard maps of the study area after each rainfall event were constructed (Figure 6). Hazard degree was classified into eight grades, and the classification was based on the method in Table 17. The I_RT of 2009 Typhoon Morakot was the highest among all rainfall events, and the number of landslide grids was 8,207; thus, the degree of landslide hazard in the study area after this typhoon was the highest. The I_RT of 2013 Typhoon Kongrey was higher than that of 0517 rainfall in 2013; thus, the degree of landslide hazard after Typhoon Kongrey was higher than that after 0517 rainfall in 2013. The I_RT of 2013 Typhoon Soulik was the lowest, and the number of landslide grids was also the lowest (n = 3,969); thus, the degree of landslide hazard after this typhoon was the lowest.

6. Conclusions

In the present study, we developed an index to evaluate landslide hazards to land use after four prominent rainfall events in Sandimen and Wutai townships in Pingtung County. Using RTC and texture analysis, we interpreted and classified the satellite images of the study area captured before and after four rainfall events. The average Kappa value was approximately 0.71, which indicated medium to high accuracy. A comparison of the satellite images captured before and after the rainfall events in terms of exposure revealed that the area of landslide due to 2009 Typhoon Morakot was the largest (1,313.12 ha); the areas of landslides due to the other rainfall events were as follows: 0517 rainfall in 2013, 813.92 ha; Typhoon Kongrey, 789.44 ha; and Typhoon Soulik, 635.04 ha. Environmental factors (e.g., elevation, slope, aspect, slope roughness, distance from the river, and geology) and slope disturbance factors (e.g., road density, building density, farmland planting rate, forest density, grassland density, and bare density) were assessed in this study to analyze slope environmental strength potential I_SESP. The average overall classification accuracy was approximately 80.4%. I_RT was calculated by multiplying EAR with I_3R,max. A new comprehensive index of rainfall-induced landslide I_RL (I_RL = I_RT × I_SESP) was established for determining the hazard of rainfall-induced landslides to land use. Using the correlation between I_SESP and I_RT, we determined the degrees of rainfall-induced landslide hazards to land use. At a constant I_SESP, higher values of I_RT indicate higher degrees of landslide hazard to land use; similarly, at a constant I_RT, higher values of I_SESP indicate higher degrees of landslide hazard to land use. Landslide occurrence is positively correlated with I_RT and I_SESP. In cases of large I_SESP values (e.g., fragile environment and high land development intensity), small I_RT values may cause landslides.