Application of High-Resolution Radar Rain Data to the Predictive Analysis of Landslide Susceptibility under Climate Change in the Laonong Watershed, Taiwan

: Extreme rainfall has caused severe road damage and landslide disasters in mountainous areas. Rainfall forecasting derived from remote sensing data has been widely adopted for disaster prevention and early warning as a trend in recent years. By integrating high-resolution radar rain data, for example, the QPESUMS (quantitative precipitation estimation and segregation using multiple sensors) system provides a great opportunity to establish the extreme climate-based landslide susceptibility model, which would be helpful in the prevention of hillslope disasters under climate change. QPESUMS was adopted to obtain spatio-temporal rainfall patterns, and further, multi-temporal landslide inventories (2003–2018) would integrate with other explanatory factors and therefore, we can establish the logistic regression method for prediction of landslide susceptibility sites in the Laonong River watershed, which was devastated by Typhoon Morakot in 2009. Simulations of landslide susceptibility under the critical rainfall (300, 600, and 900 mm) were designed to verify the model’s sensitivity. Due to the orographic e ﬀ ect, rainfall was concentrated at the low mountainous and middle elevation areas in the southern Laonong River watershed. Landslide change analysis indicates that the landslide ratio increased from 1.5% to 7.0% after Typhoon Morakot in 2009. Subsequently, the landslide ratio ﬂuctuated between 3.5% and 4.5% after 2012, which indicates that the recovery of landslide areas is still in progress. The validation results showed that the calibrated model of 2005 is preferred in the general period, with an accuracy of 78%. For extreme rainfall typhoons, the calibrated model of 2009 would perform better (72%). This study presented that the integration of multi-temporal landslide inventories in a logistic regression model is capable of predicting rainfall-triggered landslide risk under climate change. data curation, Y.-C.C. and J.-Y.T.; writing—original draft preparation, C.-W.T., C.-W.C. and C.-J.Y.; writing—review and editing, C.-W.T., C.-W.C. and C.-J.Y.; visualization, Y.-C.C.; supervision, C.-W.T. and S.-F.W.; administration, C.-W.T. and


Introduction
Taiwan is characterized by steep terrain and geological weakness co-evolving with seismic activity, monsoons, and typhoons, resulting in the high vulnerability of Taiwan. For example, the mean annual help in understanding the impacts of different types of rainfall under climate change at the landslide locations. Moreover, insights are provided for the development of follow-up adjustment strategies to improve prevention technologies and early-warning capacities for disasters in mountainous agriculture and forestry belts.

Research Area
The Laonong River is the largest tributary of the Gaoping River, with a total length of 137 km, and its watershed is the largest in Taiwan, with an area of 1409 km 2 . Most of the watershed area is within the territory of Kaohsiung City. Laonong river originates from the East peak of Yushan Mountain [26] (Figure 1). The altitude ranges from 35 to 3500 m. Geological formation includes Lushan formation, which accounts for the highest proportion of 38.9%, and the Pilushan formation, accounting for 35.1%. Therefore, shale and slate are the main geological composition. The climate is subtropical monsoon with an average annual temperature of about 24.2 °C. The annual rainfall is about 2477 to 5295 mm and precipitation of 90%, concentrated in the period from May to September [27][28][29][30][31].

Rainfall Data Sources and Classification
Rainfall data was derived from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP) and QPESUMS radar rainfall data from the Central Weather Bureau (CWB), Ministry of Transportation and Communications (MOTC). The rainfall data retrieved from TCCIP were the regular grid rainfall data with a spatial resolution of 1 × 1 km and a temporal coverage from 1960 to 2015. Rainfall data of TCCIP is produced by NCDR (National Science

Rainfall Data Sources and Classification
Rainfall data was derived from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP) and QPESUMS radar rainfall data from the Central Weather Bureau (CWB), Ministry of Transportation and Communications (MOTC). The rainfall data retrieved from TCCIP were the regular grid rainfall data with a spatial resolution of 1 × 1 km and a temporal coverage from 1960 to 2015. Rainfall data of TCCIP is produced by NCDR (National Science and Technology Center for Disaster Reduction), and MOST (Ministry of Science and Technology). NCDR collect the Remote Sens. 2020, 12, 3855 4 of 23 observation data of rainfall stations in Taiwan island scale and gridded into the 1 × 1 km resolution grid data. The weather stations include CWB automatic rainfall stations, the Water Resources Agency (WRA) MOEA (Ministry of economic affairs), and the Soil and Water Conservation Bureau (SWCB), COA (Council of Agriculture), Executive Yuan, Taiwan Power Company, and Taiwan forestry research institute, etc. The QPESUMS system is one of the remote sensing technologies and was jointly developed and built by the CWB, The NOAA (National Oceanic and Atmospheric Administration) National Severe Storms Laboratory of the United States, the Water Resources Agency (WRA) MOEA, and the SWCB, COA. The QPESUMS system integrates multiple observation data (including radar echo, rainfall station, and lightning, etc.), and combines them with a geographic information system (GIS) to enhance the real-time monitoring capability of disastrous weather such as typhoon, plum rain, afternoon convection, etc. The important development of the QPESUMS system is based on radar data. Since Doppler weather radar has the characteristics of all-weather observation and can observe large-scale rainfall and airflow characteristics, it is a very important tool for meteorological observation. We can use dense radar observations to obtain high-resolution three-dimensional (3D) data (1 × 1 km grid data) in time and space to help us understand the changes in weather systems [17]. There are 11 Doppler weather radars in Taiwan, including 4 Doppler weather radars (S-band, 10 cm wavelength, scanning every 10 min, resolution 1 km) in Hualian, Kenting, Qigu, and Zhongzheng International Airport, and 7 dual-polarization Doppler radars (C-band, 5 cm wavelength, scanning every 2 min, resolution 250 m) in Wufenshan, Linyuan, Nantun, Shulin, Qingquangang base of the air force, Magong, and Green Island. For the mountainous terrain of Taiwan, the dense radar observation network helps different radars to make up for the observation blind area caused by terrain obstruction. The related products of QPESUMS include: (1) the echo data synthesized by the radar network in Taiwan, (2) 1 and 24 h cumulative rainfall distribution-analysis of measured rainfall in the past 1 and 24 h at the whole rainfall observation station, (3) real-time monitoring of convective cells-the location of existing convective cells is detected by radar observations, (4) forecast of rain area in the next hour-based on the current rainfall situation estimated by radar data, the rain area forecast obtained by extension method, and (5) rainfall observation-rainfall values measured by rainfall stations in the past 10 min, 1, 3, 6, 12, and 24 h. Compared with the traditional method that can only use the data of ground rainfall stations to estimate the rainfall in other areas without rainfall stations, the advantage of the system is that it can objectively estimate the rainfall in areas without rainfall stations. The QPESUMS system also provides real-time radar data which can be used to accurately estimate rainfall information through comprehensive analysis [32]. The derivative products include gridded quantitative rainfall estimation and quantitative rainfall forecast products. The radar data used by the QPESUMS system in quantitative precipitation estimation is the lowest elevation angle synthesized echo that is close to the ground and is not affected by the landscape. The Z-R relationship [33] is as follows: Z = 32.5R 1.65 (1) where Z is the reflectivity and R is the rainfall rate. The estimated rainfall is corrected in cooperation with the real-time observation data from ground rain-gauge stations. Not only can the extreme values corresponding to the areas without rain-gauge stations be estimated by the system, but the estimated cumulative rainfall distribution is also more objective than the observations from a single station [34]. The available data interval is limited because of the time of system construction. The data used in this study were from 2012 to 2017. Using the QPESUMS radar rainfall data, it is possible to analyze the differences in the spatial pattern of rainfall events, as well as the rainfall hot spots.
The rainfall classification in this study is based on the criteria of heavy (torrential) rain as per the levels defined by the CWB, MOTC: heavy rain is defined as the accumulated rainfall of more than 80 mm in 24 h or hourly rainfall of more than 40 mm, extremely heavy rain is defined as an event with 24 h accumulated rainfall of more than 200 mm or 3 h accumulated rainfall of more than 100 mm, and torrential rain is defined as an event with 24 h accumulated rainfall of more than 350 mm. According to the different rainfall levels, long-term changes in the frequency of heavy rainfall in the study area were enumerated and the rainfall hot spots were analyzed.

Landslide Inventory
Satellite imagery has the advantages of a short data acquisition period, no topographic limitations, and wide image data coverage. It can be used to obtain comprehensive disaster information in real-time and is widely applied in landslide monitoring. To understand multitemporal landslides in the Laonong River watershed, the multi-temporal landslide maps and satellite images of the Laonong River from 2003 to 2015 published by the Forestry Bureau on the Taiwan Geospatial One-Stop Portal (TGOS Portal) were obtained, and we also interpreted the landslide map by using the same interpreted method [35] derived from the 2016 to 2018 SPOT (Satellite for observation of Earth) satellite imagery. Although the landslide inventory was credible and complete, the time interval of image interpretation was on the annual scale, which was too long for this study. Therefore, landslide map data were only used for analyzing the annual changes in the landslide areas, whereas the short-duration information on landslide areas after a single large-scale typhoon event was obtained from the self-interpreted SPOT satellite imagery. The normalized difference vegetation index (NDVI) was calculated by using multi-spectral information and combined with the supervised classification method, and then applied in the satellite image interpretation to explore the relationship between vegetation in landslide areas in the watershed and heavy rainfall.
NDVI has a good response to vegetation on the ground, it produces a value close to 0 in exposed soil and rock areas because they have similar reflectivity values in near-infrared and red-light bands. Because surface waters and snow-covered areas have a higher value in the red-light band than the near-infrared band, they have negative NDVI values [36]. With the above features, an NDVI value of 0.05 was used as a threshold for rapid preliminary identification in this study, with the research area preliminarily divided into vegetation and non-vegetation areas. Then, the landslide areas were extracted through manual interpretation so that the landslide area extraction could be achieved effectively and accurately. The formula for the NDVI is as follows: where NIR is the reflectivity of the near-infrared band and RED is the reflectivity of the red-light band.

Multitemporal Changes in Landslide Areas and Extraction of Explanatory Factors
Using the landslide inventory, the area of landslide sites in each year was converted to a landslide ratio, i.e., the ratio of the total area of landslides to the area of the watershed, which served as a parameter for measuring and comparing the scale of landslides in each year. In addition to rainfall, the geological and explanatory factors are also the key driving factors of landslides. The factors often considered in the analysis include elevation, slope, aspect, geology, vegetation, and land cover. With the improvement and development of instruments in recent years, high-precision elevation data can be produced through light detection and ranging (LiDAR) surveying and mapping technology. In contrast to the DEM data, the multi-penetration characteristic of LiDAR allows it to describe elevation changes under forests in mountainous areas, and the DEM data output is accurate and suitable for terrain analysis. In this study, the 20 m DEM produced by the Ministry of the Interior using LiDAR surveying was adopted to generate various factors for developing a landslide susceptibility model.

Model and Factor Selection
Multivariate logistic regression is the most commonly used method in the establishment of landslide susceptibility models [21][22][23]. Various factors including rainfall, elevation, slope, aspect, topographic curvature, topographic roughness, geological conditions, and NDVI were used as variables in the statistical regression analysis [37]. Logistic regression is a special form of log-linear model [38], and logistic distribution is a binary variable distribution function. The curve of this function is similar to the cumulative distribution curve of a random variable. When one binary variable in the log-linear model is regarded as the dependent variable and the model is defined as a series of independent variable functions, this log-linear model is a logistic regression model, whose formula is as follows: where P i = P (y i = 1|x 1i , x 2i , . . . , x ni ) is the probability of an event at the i-th grid cell when a series of independent variables (x 1i , x 2i , . . . , x ni ) are given, x ki is value of the k-th explanatory landslide factor at the i-th grid cell, and α and β k are the regression coefficients of each factor. In this study, the ten explanatory factors were used in the logistic regression. To test coefficients of the explanatory factors and their significance, we run the logistic regression model in the SPSS (Statistical Product and Service Solutions) software. The non-significant variables will be excluded from the model until all variables are significant at the 0.05 level. The calculation procedure of the ten explanatory variables is described as below: 1. Elevation (unit: meters): The altitude above sea level is directly calculated from the DEM for all of Taiwan with a 20 m resolution provided by the Ministry of the Interior.

2.
Slope (unit: degrees): Calculated from the DEM data of the Ministry of the Interior with the 20 m resolution by using ArcGIS 10.7 (Spatial Analyst of extension tool).

3.
Topographic roughness (unit: meters): The degree of topographic change in the area. Topographic roughness is calculated by using ArcGIS 10.7, with the standard deviation of the elevation values in a circular window (3 × 3 pixel, moving window) as a measure of the degree of elevation variation in the area [39]. A schematic diagram of the circular window is shown in Figure 2, and the standard deviation (SD) is calculated as follows: where SD h is the topographic roughness, n c is the number of points within the window, Z i is the elevation value of each point, and _ Z is the mean value of the elevation of each point within the window.

4.
Slope roughness (unit: degrees): The variations of slope within the area. Using the slope roughness, it is possible to find out whether a landslide would occur in the landscape with drastic slope variations. The calculation method and tool are the same as the topographic roughness, and is calculated as follows: where SD s is the slope roughness, nc is the number of points within the window, Z i is the slope value of each point, and _ Z is the mean value of the slope of each point within the window.

5.
Aspect sine value and cosine value (unitless): There will be numerical discontinuities if the aspect is expressed as an angle (e.g., 0 • = 360 • ), and thus the aspect can be decomposed into an X vector Remote Sens. 2020, 12, 3855 7 of 23 and a Y vector, i.e., the aspect sine value and aspect cosine value, respectively. A positive sine value indicates the degree to the east, and a negative sine value indicates the degree to the west, a positive cosine value indicates the degree to the north, and a negative cosine value indicates the degree to the south. The value can be calculated from the DEM data by using ArcGIS 10.7 (Spatial Analyst of extension tool). 6.
NDVI before landslides (unitless): NDVI is the most common index for evaluating the surface vegetation status in an area [40][41][42][43]. The NDVI value before the flood season can be used as a quantitative value for the vegetation status before landslides. It was calculated from SPOT satellite imagery by Erdas Imagine 4.0. 7.
Maximum accumulated rainfall in 24 h (unit: mm): In previous studies, it was often believed that rainfall intensity is related to landslide events [44,45]. However, the choice of rainfall intensity differed based on the different regions or research objectives. In previous studies on rainfall thresholds for landslides in southwestern Taiwan, the accumulated rainfall in 6, 12, 24, 48, and 72 h in autumn during different years was adopted in the variance analysis, and it was found that the 24 h accumulated rainfall had the lowest variation coefficient [46], i.e., the internal changes and differences in data corresponding to the 24 h accumulated rainfall are relatively small. The rainfall value is calculated by the Kriging method. 8.
Topographic wetness index (TWI, unit: m 2 ): The ability of the terrain to control soil moisture. The steeper the slope, the faster the runoff speed, the lower the degree of infiltration, and thereby, the lower the soil moisture content [47]. Because the watershed area is relatively large and the slope is mild in low and flat areas, the soil water content is relatively high. If the water content under the slope is large, the pore water pressure of soil increases, and the soil is affected by the lifting force of water. The value can be calculated by using ArcGIS 10.7. As a result, the probability of slope failure and sliding becomes higher, and is expressed as follows: where ω is the TWI, A s is the watershed area above a certain point, and θ is the slope. 9.
Geology (unitless): Different geological conditions affect the strength of rock mass, joint density, joint structures, etc. We used the local geological formations ( Figure 1) as a categorical variable in the model. The data are transferred from vector data to grid data. 10. Distance to the river (unit: meters): The distance from any point to the river. Landslides may occur on both sides of the riverbank. Continuous erosion on the foundation slope makes the area unstable, and a landslide is likely to occur during high-intensity rainfall. The distance from the river channel is obtained by analysis using ArcGIS 10.7 (path distance), which is to calculate the flow path from the river channel to the landslide.

Model and Factor Selection
Multivariate logistic regression is the most commonly used method in the establishment of landslide susceptibility models [21][22][23]. Various factors including rainfall, elevation, slope, aspect, topographic curvature, topographic roughness, geological conditions, and NDVI were used as variables in the statistical regression analysis [37]. Logistic regression is a special form of log-linear model [38], and logistic distribution is a binary variable distribution function. The curve of this function is similar to the cumulative distribution curve of a random variable. When one binary variable in the log-linear model is regarded as the dependent variable and the model is defined as a series of independent variable functions, this log-linear model is a logistic regression model, whose formula is as follows: where Pi = P (yi = 1|x1i, x2i, …, xni) is the probability of an event at the i-th grid cell when a series of independent variables (x1i, x2i, …, xni) are given, xki is value of the k-th explanatory landslide factor at the i-th grid cell, and α and βk are the regression coefficients of each factor. In this study, the ten explanatory factors were used in the logistic regression. To test coefficients of the explanatory factors and their significance, we run the logistic regression model in the SPSS (Statistical Product and Service Solutions) software. The non-significant variables will be excluded from the model until all variables are significant at the 0.05 level. The calculation procedure of the ten explanatory variables is described as below: 1. Elevation (unit: meters): The altitude above sea level is directly calculated from the DEM for all of Taiwan with a 20 m resolution provided by the Ministry of the Interior. Topographic roughness is calculated by using ArcGIS 10.7, with the standard deviation of the elevation values in a circular window (3 × 3 pixel, moving window) as a measure of the degree of elevation variation in the area [39]. A schematic diagram of the circular window is shown in Figure 2, and the standard deviation (SD) is calculated as follows: where ℎ is the topographic roughness, nc is the number of points within the window, Zi is the elevation value of each point, and is the mean value of the elevation of each point within the window.

Model Calibration
Logistic regression was selected as the statistical method for model construction. The size of samples in logistic regression is equal to random sampling, which is to establish the model by taking the same number of pixel samples from the landslide area and non-landslide area. The landslide map information of each year is provided by the Forestry Bureau, council of agriculture, and interpreted by satellite images. The area of non-landslide distribution in the watershed is the non-landslide area. In 2005, 5000 samples were randomly sampled in the landslide and non-landslide areas, 8000 samples were selected in the 2009 model, and 9000 samples were collected in a mixed model, including 4000 samples in 2005 and 5000 samples in 2009. Since the model is also a pixel scale analysis, sampling is also carried out at the pixel level. The environmental variables are converted into 30 m grid format data. Multiple samples may be taken from a single slope landslide site, while the environmental variables within the same landslide may be different, such as slope, aspect, or vegetative state, which affect whether the location of the pixel has high landslide susceptibility, and the identification accuracy of the calibrated model should reach over 70%.

Model Validation
The Receiver Operating Characteristic (ROC) curve is a graph used to test sensitivity and specificity, with the diagonal serving as a reference line. If the ROC curve lies on the diagonal, the landslide and non-landslide areas are not distinguishable; if the ROC curve bends toward the upper left corner from the diagonal, the test results of landslide and non-landslide areas have high sensitivity and good specificity ability [48].
Generally, in order to judge the quality of the inspection method, in addition to looking at the graph of the ROC curve, AUC (area under the ROC curve) can also be used to judge the distinguishing ability of the ROC curve. The range of the AUC value is from 0 to 1, and the larger the value means the better performance. The greater the AUC value, the higher the proportion of correct interpretation, the lower the rate of misidentification, and the better the prediction accuracy. Meanwhile, the area under the ROC curve in the graph will also be larger. Swets [49] pointed out that an AUC value of fewer than 0.5 means that the model lacks interpretability, an AUC value between 0.5 and 0.7 represents low interpretation accuracy, and a value between 0.7 and 0.9 represents high accuracy.

Landslide Risk and Critical Rainfall under Different Rainfall Scenarios
A disaster risk map is usually calculated as the product of hazard, exposure, and vulnerability [50]. In the Laonong watershed, the land-uses are primarily indigenous forests and rarely cultivated fields and building areas. The importance of the forest resource and human-made facilities may be different for different management authorities (e.g., disaster prevention and forest management authorities). Besides, the occurrence of the landslide can destroy all human life and infrastructures, so the vulnerability of the facilities can be ignored. Therefore, the factors of exposure and vulnerability were ignored, and only the factor of the hazard was considered in this risk map. It should be noted that the simplification of the risk factors may underestimate the landslide risk in the building areas for disaster management. For the application of the best model, we assume that the geological and topographical factors are fixed, and the NDVI can only be generated using the latest (2019) image. Therefore, rainfall is used as the designated factor to simulate landslide susceptibility under the different rainfall intensities. The probability of landslide under different rainfall scenarios was estimated based on the explanatory parameters, and the simulated probability values were classified into five categories (Table 1). In addition, the landslide susceptibility models were used to infer critical rainfall that induces disasters, and the critical rainfall values in the study area under the disaster probabilities of 50%, 70%, and 90% were simulated, which can be used as standards for rainfall warnings in the future to prevent disasters before they happen.

Long-Term Rainfall Variations
According to the long-term changes in annual rainfall of TCCIP from 1960 to 2017 (Figure 3a (Figure 3c), it can be seen that the frequency after 2004 was significantly higher than that in the previous 40 years, which may be related to the number of typhoons. After 2000, the path tracking of typhoons in the northwest Pacific Ocean generally shifted northward, causing significantly more typhoons to invade Taiwan [51]. Moreover, the probability of extremely heavy rainfall events related to typhoons has also increased greatly [52]. Compared with the increase of annual rainfall after 2004 and the frequency change of heavy rain, it shows that the increase of annual rainfall is closely related to the occurrence frequency of heavy rain, that is, the main contribution of the increased annual rainfall comes from the form of heavy rainfall. Besides, the TCCIP and QPE sums exhibit significant differences in annual rainfall (Figure 3a) but are corresponding in the events of daily rainfall up to 80 and 200 mm (Figure 3b,c). This underestimation of TCCIP in annual rainfall resulted from the relatively few rainfall stations and missed some local precipitation in the mountainous areas.

The Multitemporal Changes in Landslide Areas
The annual changes in the area of landslides in the Laonong River watershed from 2003 to 2018 are shown in Table 2. From 2003 to 2008, the area of landslides was less than 2700 ha, and the landslide ratio of the watershed was about 1.5%. The statistics show that after Typhoon Morakot (2009), the area of landslides increased to 9472.67 ha, and the landslide ratio of the watershed rose to 6.96%. Compared with the average landslide ratio of 1.5% before 2008, approximately 5.46% of land cover was transformed into landslide areas.

The Multitemporal Changes in Landslide Areas
The annual changes in the area of landslides in the Laonong River watershed from 2003 to 2018 are shown in Table 2. From 2003 to 2008, the area of landslides was less than 2700 ha, and the landslide ratio of the watershed was about 1.5%. The statistics show that after Typhoon Morakot (2009), the area of landslides increased to 9472.67 ha, and the landslide ratio of the watershed rose to 6.96%. Compared with the average landslide ratio of 1.5% before 2008, approximately 5.46% of land cover  After Typhoon Morakot, the area of landslides showed a gradually declining trend ( Figure 4). Specifically, the area of landslides reduced sharply from 2009 to 2012, indicating that the landslide areas caused by Typhoon Morakot had recovered, but there was still a long way to go to reach the status characterizing the area before the typhoon disaster. The area of landslides from 2013 to 2018 remained between 4786 and 6352 ha, and the landslide ratio of the watershed fluctuated between 3.5% and 4.5%. The change in the landslide area ratio shows that, after Typhoon Morakot in 2009, vegetation cover was restored in a considerable amount of landslide area within 3 years. However, the recovery from large-scale deep-seated landslides is difficult. During the rainy season, repeated landslides prevented the landslide ratio from reaching the level before Typhoon Morakot in 2009.

Landslide Susceptibility Model and Explanatory Factors
We first used the landslide inventory of 2009 Typhoon Morakot to establish the landslide susceptibility model. The results of logistic regression analysis (Table 3) show that there are nine optimal explanatory factors: the maximum continuous rainfall in 24 h, elevation, slope, aspect sine value, aspect cosine value, TWI, distance to the river, NDVI before landslides, and geology. Moreover, the overall accuracy from the prediction model reached 72.12% (Table 4), and the AUC value of 0.788. Then, the tested model of 2009 is further verified by using the landslide of 2005 Typhoon Haitang. However, the AUC value declined to 0.639 (Table 4). This result indicates that a model calibrated using an intense rainfall event may not be capable of predicting the landslide triggered by moderate rainfall. Before Typhoon Morakot in 2009, the average landslide area was less than 1.28 ha, which increased to 2.05 ha after the typhoon, showing that many large-scale landslides occurred. From 2009 to 2010, the number of landslide sites increased from 4615 to 7337, and the landslide area decreased by 1005 ha. In 2010, the average landslide area was reduced to 1.09 ha, i.e., the number of landslide sites increased whereas the area of landslides decreased. One possible reason could be that the invasion and growth of herbaceous plants on areas affected by large-scale landslides aided in the recovery of the local vegetation and fragmented the shape of single large landslides. Thus, the site of a previously large landslide can be counted as multiple landslide areas in the interpretation. The average landslide area between 2011 and 2018 remained at 0.85-1.67 ha, except for the slight increase in the landslide scale due to the successive impacts of Typhoon Nibert and Typhoon Megi in 2016. In general, after Typhoon Morakot, the average landslide area in the Laonong River watershed remained higher than that before the typhoon disaster for a long time, which might be related to the fact that land recovery is difficult in large-scale landslide areas.

Landslide Susceptibility Model and Explanatory Factors
We first used the landslide inventory of 2009 Typhoon Morakot to establish the landslide susceptibility model. The results of logistic regression analysis (Table 3) show that there are nine optimal explanatory factors: the maximum continuous rainfall in 24 h, elevation, slope, aspect sine value, aspect cosine value, TWI, distance to the river, NDVI before landslides, and geology. Moreover, the overall accuracy from the prediction model reached 72.12% (Table 4), and the AUC value of 0.788. Then, the tested model of 2009 is further verified by using the landslide of 2005 Typhoon Haitang. However, the AUC value declined to 0.639 (Table 4). This result indicates that a model calibrated using an intense rainfall event may not be capable of predicting the landslide triggered by moderate rainfall.  Due to the low accuracy in the calibration model in 2009, the study also used the 2005 Typhoon Haitang for developing the model. The optimal factor combination of this model included eight factors: the maximum continuous rainfall in 24 h, slope, aspect sine value, aspect cosine value, TWI, distance to the river, NDVI before landslides, and geological conditions, and the overall accuracy from the prediction model could reach 77.78% and AUC values of 0.895. By using the landslide of 2009 Typhoon Morakot, the AUC value, however, declined to 0.619. The 2005 calibrated model applies not only to rare and extreme events but also to most typhoon-induced landslide events.
This study aimed to establish a general landslide susceptibility model. Therefore, a mixed sampling of the two typhoon-induced landslide events was adopted for model establishment. The overall prediction accuracy of the calibrated model was 69.81% (Table 4), and the optimal factor combination included ten factors: the maximum continuous rainfall in 24 h, elevation, slope, topographic roughness, aspect sine value, aspect cosine value, TWI, distance to the river, NDVI before landslides, and geology. The performance was not as good as that of the 2005 calibrated model. However, its validation performances on both events were acceptable (0.772 for 2005 Typhoon Haitang and 0.811 for 2009 Typhoon Morakot), showing that the mixed calibrated model described in this study can be used as an auxiliary tool for the general model.
By comparing the Wald value performances of the three sets of calibrated models, the controlling factors in each model could be determined. In the 2005 calibrated model, the NDVI before landslides and geology are the main controlling factors. The NDVI before landslides plays a critical role, and its odds ratio (Exp β) is 0.001, indicating that vegetation coverage helps to reduce the probability of landslides under the background of perennial typhoon-induced landslide events. In the 2009 calibrated model, the main controlling factors are the aspect of cosine value, aspect sine value, geology, and TWI. Among these, the aspect factors are the key parameters. Although the NDVI before landslides reached a significant level, its Wald value was the lowest among the group of parameters, implying that under large-scale extreme rainfall events, the impact of vegetation coverage on slope stability is significantly reduced and that the topographic conditions become dominant. The important parameters in the mixed calibrated model are the aspect cosine value, NDVI before landslides, aspect sine value, geology, and TWI, including key parameter adjustments of vegetation and topographic conditions. As a result, the model validation performance during extreme events and perennial typhoon events is not as good as those of the other single-event calibrated models.

Rainfall Hot Spots and Landslides
The changes in the long-term rainfall frequency and the spatial distribution of rainfall in the research area were analyzed based on the QPESUMS data from 2012 to 2017 ( Figure 5). It can be seen that the rainfall is concentrated in two hot spots. The southern rainfall center lies at the junction of the Liouguei District and Maolin District in Kaohsiung City and extends northward to the southern part of the Taoyuan District, while the northern rainfall center is in the northern part of Taoyuan District, located in the upstream valley of the Laonong River. The southern rainfall center is larger than the northern center in both the area coverage and the rainfall amount.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 24 landslides under the background of perennial typhoon-induced landslide events. In the 2009 calibrated model, the main controlling factors are the aspect of cosine value, aspect sine value, geology, and TWI. Among these, the aspect factors are the key parameters. Although the NDVI before landslides reached a significant level, its Wald value was the lowest among the group of parameters, implying that under large-scale extreme rainfall events, the impact of vegetation coverage on slope stability is significantly reduced and that the topographic conditions become dominant. The important parameters in the mixed calibrated model are the aspect cosine value, NDVI before landslides, aspect sine value, geology, and TWI, including key parameter adjustments of vegetation and topographic conditions. As a result, the model validation performance during extreme events and perennial typhoon events is not as good as those of the other single-event calibrated models.

Rainfall Hot Spots and Landslides
The changes in the long-term rainfall frequency and the spatial distribution of rainfall in the research area were analyzed based on the QPESUMS data from 2012 to 2017 ( Figure 5). It can be seen that the rainfall is concentrated in two hot spots. The southern rainfall center lies at the junction of the Liouguei District and Maolin District in Kaohsiung City and extends northward to the southern part of the Taoyuan District, while the northern rainfall center is in the northern part of Taoyuan District, located in the upstream valley of the Laonong River. The southern rainfall center is larger than the northern center in both the area coverage and the rainfall amount. The main reasons for that are the influence of concentrated rainfall in summer and plum rain season, and the influence of terrain amplitude effect. In summer, typhoon and southwest monsoon The main reasons for that are the influence of concentrated rainfall in summer and plum rain season, and the influence of terrain amplitude effect. In summer, typhoon and southwest monsoon bring strong wind and rainfall. In the shallow mountain area and middle altitude area in the south of Laonong River watershed, a large amount of rainfall is caused by the combination of terrain amplitude. The rainfall decreases in the middle part of the river valley, and some water vapor can be intercepted near the high altitude of the central mountain range. In autumn, the rainfall in the watershed decreased significantly, and the rainfall gradually shifted to the high-altitude area, which also showed that the seasonal rainfall variation in the southwest of the watershed was quite large. For example, southwest-facing slopes are the windward side during typhoon Morakot, and typhoon-induced landslides are concentrated on southeast-facing slopes [53]. In contrast, the southeast slopes are mostly forward slopes. Heavy rainfall penetrates between the rock formation and increases the pore water pressure, which may increase the probability of landslide by reduced friction [54].
To understand the spatial distribution of extremely heavy rainfall, the top 1% of rainfall events were selected based on the hourly rainfall and the 3 h accumulated rainfall to map the spatial distribution of extreme rainfall intensity. The results show that the spatial distribution patterns of heavy rain and extremely heavy rainfall events are the same. The top 1% of the average rainfall is mainly distributed in a hot spot in the south of the watershed. The hourly rainfall intensity of these events can generally reach more than 30 mm (Figure 6a), and the accumulated rainfall in 3 h can reach 65 mm (Figure 6b). This hot spot can be further divided into two rainfall centers, the valley area in the low mountains on the southwest side, and the central-southern watershed, along the dividing crest of the Baolai.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 24 bring strong wind and rainfall. In the shallow mountain area and middle altitude area in the south of Laonong River watershed, a large amount of rainfall is caused by the combination of terrain amplitude. The rainfall decreases in the middle part of the river valley, and some water vapor can be intercepted near the high altitude of the central mountain range. In autumn, the rainfall in the watershed decreased significantly, and the rainfall gradually shifted to the high-altitude area, which also showed that the seasonal rainfall variation in the southwest of the watershed was quite large. For example, southwest-facing slopes are the windward side during typhoon Morakot, and typhooninduced landslides are concentrated on southeast-facing slopes [53]. In contrast, the southeast slopes are mostly forward slopes. Heavy rainfall penetrates between the rock formation and increases the pore water pressure, which may increase the probability of landslide by reduced friction [54].
To understand the spatial distribution of extremely heavy rainfall, the top 1% of rainfall events were selected based on the hourly rainfall and the 3 h accumulated rainfall to map the spatial distribution of extreme rainfall intensity. The results show that the spatial distribution patterns of heavy rain and extremely heavy rainfall events are the same. The top 1% of the average rainfall is mainly distributed in a hot spot in the south of the watershed. The hourly rainfall intensity of these events can generally reach more than 30 mm (Figure 6a), and the accumulated rainfall in 3 h can reach 65 mm (Figure 6b). This hot spot can be further divided into two rainfall centers, the valley area in the low mountains on the southwest side, and the central-southern watershed, along the dividing crest of the Baolai.    (Table 5) show that under the rainfall scenario of 300 mm, the area characterized at the high-risk level or above accounted for 3.95% of the watershed, which gradually increased to 5.07% under 600 mm rainfall conditions and 6.68% under 900 mm rainfall. The simulation results of the mixed model (Table 6) show that under the rainfall scenario of 300 mm, the area corresponding to the high-risk level or above accounted for 4.51% of the watershed, which increased to 5.20% under 600 mm rainfall and 6.02% under 900 mm rainfall. In 2005, the maximum 24 h average rainfall in the watershed caused by Typhoon Haitang was approximately 667 mm, and the landslide ratio of the watershed was 1.96%. Comparing the simulation results for areas of high-risk level or above under a rainfall scenario of 600 mm given by the two models, the landslide area was overestimated by approximately 3%. For Typhoon Morakot in 2009, the maximum 24 h average rainfall was approximately 1020 mm, and the landslide ratio of the watershed after the disaster was 6.96%. Comparing the simulation results for areas of high-risk level or above under a rainfall scenario of 900 mm given by the two models, the landslide area was slightly underestimated by 0.28% (2005 model) and 0.94% (mixed model). Although there were slight differences observed in the comparison between the landslide ratios of historical disaster events and the 2019 parameter simulations, which are mainly comparison errors due to the different conditions of vegetation factors, the models can be used to measure the actual differences in landslides induced by historical disastrous rainfall events. The simulations with the two models revealed that most areas in the watershed are characterized by slight landslide risk. However, the mixed calibrated model gives a higher proportion of areas with low-to-moderate risk levels and the proportion of areas with low risk or above increases when rainfall is lower than that in the simulation results of the 2005 calibrated model, which is related to the contribution of topographic conditions in the mixed model. In the application of risk ranking, if a single event is estimated to be high-intensity homogeneous rainfall, the mixed model provides better results. In this study, however, our goal is to present that models calibrated from landslide inventories with different rainfall intensity have different weight of explanatory variable. For the test model in 2005, the previous NDVI is an important explanatory variable, but not for the test model in 2009 (Table 3). This could be because great intensive rainfall in 2009 Typhoon Morakot decreased the shear strength of hillslopes which have already largely exceeded the effects of root cohesion reinforcement. Therefore, the results suggest that a mixed calibrated model could be used as a supplement to the general model but cannot easily catch all the landslide with different rainfall intensity. For example, Chang et al. suggested that a catch-all model (combine all landslide datasets with different rainfall intensity) can predict landslides triggered by a strong typhoon (intensive rainfall intensity) but not a weak typhoon (low rainfall intensity), because a weak typhoon usually triggers a small number of landslides in isolated areas, which can be missed by a model developed from other typhoon events [55]. With fixed topographic conditions, the possible critical values under different disaster probabilities can be calculated using the logistic regression model. Three levels of landslide probability, i.e., 50%, 70%, and 90%, were used in the simulations to analyze the accumulated rainfall in 24 h that may cause disasters in each area. The study results can be used as a reference in the form of regional disaster-prevention maps for villages and residents (Figures 7 and 8). During rainstorm and typhoon events, the evaluation results based on the regional observation values can be applied to prevent disasters from occurring.

Conclusions
Spatial distribution of annual rainfall in the Laonong River watershed indicated two rainfall hot spots, the junction of the Liouguei District and Maolin District in Kaohsiung City, and the northern part of Taoyuan District. Additionally, frequency of heavy rain and extremely heavy rain has increased significantly, which may suggest that satellite images should be updated before the flood season every year for simulating the disaster risk. According to the changes in landslide areas, the average landslide ratio before 2008 was 1.5%, which rapidly increased to 6.96% after Typhoon Morakot and fluctuated between 3.5% and 4.5% after 2013. The average landslide area remained at a

Conclusions
Spatial distribution of annual rainfall in the Laonong River watershed indicated two rainfall hot spots, the junction of the Liouguei District and Maolin District in Kaohsiung City, and the northern part of Taoyuan District. Additionally, frequency of heavy rain and extremely heavy rain has increased significantly, which may suggest that satellite images should be updated before the flood season every year for simulating the disaster risk. According to the changes in landslide areas, the average landslide ratio before 2008 was 1.5%, which rapidly increased to 6.96% after Typhoon Morakot and fluctuated between 3.5% and 4.5% after 2013. The average landslide area remained at a high level for a long time, which is speculated to be related to the difficulty in restoration after large-scale landslides.
Calibration and validation results of the landslide susceptibility models showed that the calibrated model of 2005 is preferred for conditions of perennial typhoon-induced landslides. For extremely heavy rainfall, the calibrated model of 2009 should be adopted for better simulation results. The results of the logistic regression model showed that the vegetation before rainfall event is related to the reduction of landslide susceptibility. It indicates that the treatment of landslide in the Laonong River watershed is effective. The assessment of the landslide susceptibility model in this study can draw up a long-term landslide-control strategy for rainstorm events, which can improve the efficiency of landslide area remediation. For extreme rainfall events with a high risk of disaster, disaster damage assessment can also be carried out in advance, and action strategies can be planned to achieve the purpose of disaster reduction.