A Spatial Model of Landslides with a Micro-Topography and Vegetation Approach for Sustainable Land Management in the Volcanic Area

: This study aims to produce a spatial model for sustainable land management in landslide-prone areas, based on exploring non-stationary relationships between landslide events, geomorpho-logical and anthropogenic variables on tropical hillsides, especially in Taji Village, Jabung District, East Java Province, Indonesia. A series of approaches combine in this research, and methods are used to construct independent and dependent variables so that GWR can analyze them to obtain the best model. Transformation of categorical data on microtopography, landform, and land cover variables was carried out. When modelled, landscape metrics can explain landslide events in the study area better than distance metrics with adj. R 2 = 0.75 and AICc = 2526.38. Generally, local coefﬁcient maps for each variable are mapped individually to reveal their relationship with landslide events, but in this study they are integrated to make it more intuitive and less confusing. From this map, it was found that most of the variables that showed the most positive relationship to the occurrence of landslides in the study area were the divergent footslopes. At the same time, the negative one was plantation land. It was concluded that the methodological approach offered and implemented in this study provides signiﬁcant output results for the spatial analysis of the interaction of landslide events with geomorphological and anthropogenic variables locally, which cannot be explained in a global regression. This study produces a detailed scale landslide-prone conservation model in tropical hill areas and can be reproduced under the same geo-environmental conditions.


Introduction
Land degradation for several countries is one of the problems that can lead to disasters [1,2] by causing a reduction or loss of land productivity, resulting in economic losses. Based on [3] 25% of land area worldwide is degraded, and as much as 24 billion tons of fertile soil are lost every year due to degradation. Land degradation not only causes disasters, but disasters also cause land degradation, one of which is landslides.
Landslides are natural disasters that usually occur in mountainous or hilly areas. These disasters often cause extensive economic losses and yearly fatalities [4]. Indonesia is located above the confluence of three major plates, namely the Eurasian plate, the Pacific plate and the Indo-Australian plate, reflecting high tectonic activity with a tropical climate and intensive anthropogenic activity, which often causes natural disasters. According to the Indonesian Disaster Information Data, landslides are ranked 3rd (9047 incidents) after tornadoes (11,016 incidents) and floods (13,723 incidents) recorded since 1822 until now [5]. Landslide disasters in Java Island, Indonesia tend to be caused by high rainfall, which often occurs in remote hilly areas that are prone to these events [6,7]. In the vicinity of Mount Bromo, landslides occur due to precipitated volcanic material, steep slopes, and high rainfall which often damage road accessibility [8]. Landslides caused by high rainfall are a global problem yearly [9]. Landslides in Taji Village-one of the villages in the Mount Bromo area-are caused by high-intensity rains and extreme weather (La nina), and often damage the connecting roads between villages which causes residents to be isolated, thus hindering farmer and gardening activities in the fields. In addition, the houses of affected residents have also occurred, but there were no fatalities.
In addition to high rainfall as a trigger, the phenomenon of slope instability that causes landslides can also be affected by landform conditions in terms of morphology, morphoprocess, morphochronology, and morphoarrangement [10]. In addition, geomorphological mapping in mountainous and hilly areas is the most complex type of information and the most subjective for landslide hazard assessment [11]. However, geomorphological mapping is selective because it only focuses on certain features in an area of interest with a certain scale for a particular study. The latest remote sensing techniques, namely by utilizing Unmanned Aerial Vehicles (UAV) data and Digital Terrain Models (DTM) can expand their application in geomorphological and topographical mapping for the identification and mapping of landslide hazards [12][13][14] in more detailed area coverage at a precise scale. In addition, recently, low-cost UAVs (more commonly called drones) have become a trend among academics, practitioners and commercial circles because they are effective in collecting large amounts of elevation data in a relatively short time, which can change the perspective and analysis by geomorphologists to study geomorphometrics and topography in certain landscapes [15,16].
One of the geomorphological features on a detailed scale is the microtopography built from the DTM. Micro-topography is defined as topographic changes on a small scale that is divided into seven units, including crest slope, upper side slope, head hollow, lower side slope, flood terrace and riverbed [17,18]. Chimner dan Hart (1996) defines microtopography into three units: hummock, pool and intermediate area. In terms of the scope of soil development, microtopography is divided into two types: pit and mound on soil morphology caused by fallen trees forming drumlin landforms [19,20]. Another definition of microtopography is the difference in size and shape in the local terrain caused by soil erosion, thus affecting the heterogeneity of habitat conditions such as moisture and soil nutrients on a scale of 1 m 2 or more [21,22] with Microtopographic types include platform, gully, sink hole, scarp and ephemeral gully [23]. Thus, the definition and type of microtopography is "variable", which adjusts to the study in a particular field. In this study, microtopography is defined as changes in local topography in terms of size and shape, as seen from differences in morphology [24,25]. Morphological mapping is based on line shape mapping, which focuses on identifying the types of slope bends using a symbology system that is unambiguous, clear, and reproducible [26]. In other words, microtopography is a reflection of its morphology.
In general, geomorphological features have a major influence on landslide events. Anthropogenic activity also plays an important role in slope instability, part of which is land cover [27,28], and even contributes positively to landslide events [29]. However, land cover and geomorphological features are closely related, allowing different spatial relationships to landslide events [30]. In addition, vegetation density can also explain the pattern of landslide occurrence, but often has an inverse relationship, namely the higher the vegetation density, the lower the landslide vulnerability [13,31].
Spatial modeling of landslide susceptibility is crucial for further understanding the assessment of disaster mitigation and preventive measures to conserve land in landslideprone areas. Approaches for landslide hazard mapping are developing rapidly, starting from combined methods to produce the best models to emerging innovations for model updates. In general, they include the evaluation of landslide models using qualitative approaches [32], quantitatively based on the relationship of controlling factors and landslides [33], and even a combination of both [34,35]. However, the spatial interaction between landslide points and their controlling factors in the quantitative approach is not explained [35], it only relies on stationary parameter estimation to examine the relationship between the two. One of the quantitative approaches, local regression analysis, can explain these spatial interactions in a non-stationary manner [36]. Feuillet et al.'s (2014) study examined the strength of the spatial relationship between paraglacial factors and landslide events. In addition, the authors of [29] investigated the local spatial relationship between land use change resulting from human intervention and landslide events. Research [24] also explains that different vegetation classes on different microtopography give different responses to the process of soil loss in the form of erosion. These studies prove that the occurrence of landslides has a non-stationary relationship to the predictor variable, which is at the same time better than global regression in general. However, the coverage area of the two studies is on a regional scale.
In this context, there are rarely studies that discuss the spatial relationship between landslide events and landslide control factors at a detailed scale, especially in the tropical hills of Indonesia. In fact, the spatial relationship of landslides with their controlling factors can provide information on the biggest factors causing landslides spatially. Consequently, it will be known that the arrangement of the microtography and vegetation parts that have a positive and negative effect on landslides. This information can be essential information in the management of sustainable land management.
As previously explained, drones are currently becoming a trend and are applied for specific purposes, as well as an alternative to optical satellite imagery data with very high spatial resolution that are quite expensive, and user demands are also high. Thus, this study aims to produce a detailed scale conservation model in landslide-prone areas based on exploring the local spatial relationship between landslide events and micro-topographical variables, land cover, and vegetation density at a detailed scale in a small hilly area in Taji Village, Jabung District, Province East Java, Indonesia uses the Geographically Weighted Regression (GWR) method.

Study Area
The study area is in the Bromo Tengger Semeru Area ( Figure 1). Astronomically, it is located at 7 • 56 34.98'-7 • 57'6.1" South Latitude and 112 • 48'49"-112 • 49'30.58" East Longitude. The study area covers 61.2 ha with an average elevation of 1110 ± 59 masl. The topographical characteristics in the study area are hillsides that are quite steep to steep and cut by rivers to form a fairly deep valley. The use of agricultural land and plantations tends to dominate in the study area. On land that is not vegetated erosion is found to be more intensive. The geological conditions in the study area consist of lower quarter volcanic rocks (i.e., Mount Gendis) during the middle Pleistocene. The rock materials include volcanic breccias, tuff-breccias, lava, and agglomerates. In addition, areas with andesitic rock deposits, namely lava and breccia-tuff, tend to be prone to landslides [37]. The complex topography configuration, high erosion rate, and rock materials in the study area have great potential for future landslides.

Methodology
Exploration of spatial relationships locally using the GWR model between landslides and microtopographical variables, landform, land cover, and vegetation density through several stages, involved: (1) building a Digital Terrain Model (DTM) from overlapping drone photos in the study area; (2) creating imagery orthomosaic based on dense point cloud, mesh, and texture data; (3) preparing raw data in the form of orthomosaic imagery acquired from UAV drones and landslide inventory via orthomosaic; (4) orthomosaic imagery used for land cover analysis using the Geographically Object-based Image Analysis (GEOBIA) method) based on spectral features, haralick texture, and shape index, which then attribute selection is carried out for all features through WEKA software to produce optimal land cover classification; (5) individual stand identification from orthomosaic image interpretation for vegetation density analysis; (6) curvature-based landform classification by classification system Pennock uses DTM; (7) micro-topography constructed through on-screen digitization based on elevation contour lines and the landscape appearance of the study area, which is then zoning micro-topography using the Voronoi diagram or Thiessen polygon method; (8) integrating landslide data, landform, land cover, and vegetation densityinto microtopography zoning as a spatial unit based on the value of the results of transforming categorical data into numeric (specifically for microtopography, landform, and land cover); (9) using three types of model for each variable, namely Type I model (proximity factor), Type II (Principal Component (PC) with the highest percentage of eigenvalues on the landscape metric comprehensive index), and Type III (PC on the landscape metric comprehensive index with the largest contribution using the Relief-F attribute selection method from WEKA software) prepared for GWR model fitting; and 10) local spatial relationship analysis based on the best fit GWR model with the four variables bell is modeled simultaneously. This series of stages can be simplified through the research flowchart in (Figure 2).

Methodology
Exploration of spatial relationships locally using the GWR model between landslides and microtopographical variables, landform, land cover, and vegetation density through several stages, involved: (1) building a Digital Terrain Model (DTM) from overlapping drone photos in the study area; (2) creating imagery orthomosaic based on dense point cloud, mesh, and texture data; (3) preparing raw data in the form of orthomosaic imagery acquired from UAV drones and landslide inventory via orthomosaic; (4) orthomosaic imagery used for land cover analysis using the Geographically Object-based Image Analysis (GEOBIA) method) based on spectral features, haralick texture, and shape index, which then attribute selection is carried out for all features through WEKA software to produce optimal land cover classification; (5) individual stand identification from orthomosaic image interpretation for vegetation density analysis; (6) curvature-based landform classification by classification system Pennock uses DTM; (7) micro-topography constructed through on-screen digitization based on elevation contour lines and the landscape appearance of the study area, which is then zoning micro-topography using the Voronoi diagram or Thiessen polygon method; (8) integrating landslide data, landform, land cover, and vegetation densityinto microtopography zoning as a spatial unit based on the value of the results of transforming categorical data into numeric (specifically for microtopography, landform, and land cover); (9) using three types of model for each variable, namely Type I model (proximity factor), Type II (Principal Component (PC) with the highest percentage  Very high-resolution aerial photos were acquired through photogrammetric processing using a UAV (a multirotor-type drone, the DJI Phantom 4). The Pix4Dcapture software [38] is used for automatic flight control and aerial photo acquisition to retrieve information on surface objects in the study area. In addition, geometric correction of orthorectified images is also automatically performed in Pix4D. The mapping was carried out at an altitude of about 70 m, producing an image with a spatial resolution of 2.4 cm per pixel. For the flight path, forward and side overlap when shooting is set optimally at around 80% and 70%, respectively. This is due to reducing the canopy height error in vegetation, which is a larger proportion than non-vegetation in the study area [39]. The aerial photos that have been acquired are then processed using third-party software Agisoft PhotoScan [40] to build orthomosaic images and Digital Terrain Models (DTM), which are often used in previous research for experimental and other scientific fields, especially land cover mapping [39,[41][42][43][44][45].  Very high-resolution aerial photos were acquired through photogrammetric processing using a UAV (a multirotor-type drone, the DJI Phantom 4). The Pix4Dcapture software [38] is used for automatic flight control and aerial photo acquisition to retrieve information on surface objects in the study area. In addition, geometric correction of orthorectified images is also automatically performed in Pix4D. The mapping was carried out at an altitude of about 70 m, producing an image with a spatial resolution of 2.4 cm per pixel. For the flight path, forward and side overlap when shooting is set optimally at around 80% and 70%, respectively. This is due to reducing the canopy height error in vegetation, which is a larger proportion than non-vegetation in the study area [39]. The aerial photos that have been acquired are then processed using third-party software Agisoft PhotoScan [40] to build orthomosaic images and Digital Terrain Models (DTM), which are often used in previous research for experimental and other scientific fields, especially land cover mapping [39,[41][42][43][44][45]].

Landslide Inventory
The landslide inventory map was compiled from orthomosaic imagery, because historical data on past landslide events were not recorded. High-resolution imagery from Google Earth with a timestamp also cannot help to identify landslide events because the scope and scale in the study area is very detailed. A total of 14 landslide points were detected in the study area ( Figure 1). Then these landslide data are aggregated into microtopographic units as the dependent variable for the GWR model. Generally, landslides are denoted as a binary class, namely, 0 (not landslide) and 1 (landslide). However, this does not represent the actual landslide events when overlaid onto microtopographic units. In other words, overlapping landslide areas do not always completely intersect with microtopographic units because these spatial units certainly form non-uniform areas such as grids with a fixed area shape. In addition, a statistical model for slope instability can use the percentage of landslide area in each unit of analysis [33].

Landslide Inventory
The landslide inventory map was compiled from orthomosaic imagery, because historical data on past landslide events were not recorded. High-resolution imagery from Google Earth with a timestamp also cannot help to identify landslide events because the scope and scale in the study area is very detailed. A total of 14 landslide points were detected in the study area ( Figure 1). Then these landslide data are aggregated into micro-topographic units as the dependent variable for the GWR model. Generally, landslides are denoted as a binary class, namely, 0 (not landslide) and 1 (landslide). However, this does not represent the actual landslide events when overlaid onto microtopographic units. In other words, overlapping landslide areas do not always completely intersect with microtopographic units because these spatial units certainly form non-uniform areas such as grids with a fixed area shape. In addition, a statistical model for slope instability can use the percentage of landslide area in each unit of analysis [33].

Microtopographic Zoning
The morphology to be mapped refers to the basic classification of Cooke dan Doornkamp (1974), including cliff, an angular convex break of slope, an angular concave break of slope, smoothly convex change in slope, smoothly concave change in slope, and convex and concave too close together (breaks of slope and smooth change in slope, respectively). However, the use of this classification is slightly subjective as there are no definitive rules as to what an angular or smooth break of slope actually looks like, the extent to which small undulations (<1 m) can or should be mapped, and where a break of slope no longer occurs [26]. Apart from that, this research may change or add morphological classes based on the landscape characteristics in the study area. The morphological mapping process utilizes elevation contour lines from the DTM data that have been made. The on-screen image interpretation technique [10,46] is used for morphological delineation through observing contour line patterns and orthomosaic imagery to see the characteristics of the landscape. Additionally, the output from the results of morphological mapping is used for micro-topographical zoning.
However, the spatial form of morphology is in the form of vector lines. Microtopographic zoning should be in the form of areas or polygons so that an approach is needed to convert the form of spatial data. The buffer zone approach defines the area's boundaries within the morphological unit to the adjacent morphology. In fact, each unit's buffer zones can overlap because the buffer distance is fixed [47]. The expected output is a flexible buffer zone that does not coincide, meaning that microtopographic zoning will be formed when the buffer zone boundary in the morphological unit touches the buffer zone boundary of the adjacent unit. Therefore, the Tyson polygon technique, also called the Thiessen polygon, is used to overcome the problem of overlapping buffer zones. The morphological output results are converted into points because vector lines consist of more than one vertex point, so the Thiessen polygon technique can be executed. Then, each point vertex whose area has been formed is aggregated based on the same ID, namely the morphological unit. For more details, see the schematic diagram for microtopographic zoning in (Figure 3).
small undulations (<1 m) can or should be mapped, and where a break of slope no longer occurs [26]. Apart from that, this research may change or add morphological classes based on the landscape characteristics in the study area. The morphological mapping process utilizes elevation contour lines from the DTM data that have been made. The on-screen image interpretation technique [10,46] is used for morphological delineation through observing contour line patterns and orthomosaic imagery to see the characteristics of the landscape. Additionally, the output from the results of morphological mapping is used for micro-topographical zoning.
However, the spatial form of morphology is in the form of vector lines. Micro-topographic zoning should be in the form of areas or polygons so that an approach is needed to convert the form of spatial data. The buffer zone approach defines the area's boundaries within the morphological unit to the adjacent morphology. In fact, each unit's buffer zones can overlap because the buffer distance is fixed [47]. The expected output is a flexible buffer zone that does not coincide, meaning that microtopographic zoning will be formed when the buffer zone boundary in the morphological unit touches the buffer zone boundary of the adjacent unit. Therefore, the Tyson polygon technique, also called the Thiessen polygon, is used to overcome the problem of overlapping buffer zones. The morphological output results are converted into points because vector lines consist of more than one vertex point, so the Thiessen polygon technique can be executed. Then, each point vertex whose area has been formed is aggregated based on the same ID, namely the morphological unit. For more details, see the schematic diagram for microtopographic zoning in (Figure 3).

Landforms of the Pennock Classification System
The landform analyzed in this study includes the landform elements based on the surface shape on the slopes of the hills, which are explained by topographical derivatives, namely slope, plan curvature and profile curvature proposed by [48]. DTM data are used for the classification of landforms that include: Convergent Shoulder (LF1); Convergent Backslope (LF2); Convergent Footslope (LF3); Divergent Shoulder (LF4); Divergent Backslope (LF5) and Divergent Footslope (LF6). Each of these landform classes is identified based on the threshold value of the combination of degrees of slope and curvature

Landforms of the Pennock Classification System
The landform analyzed in this study includes the landform elements based on the surface shape on the slopes of the hills, which are explained by topographical derivatives, namely slope, plan curvature and profile curvature proposed by [48]. DTM data are used for the classification of landforms that include: Convergent Shoulder (LF1); Convergent Backslope (LF2); Convergent Footslope (LF3); Divergent Shoulder (LF4); Divergent Backslope (LF5) and Divergent Footslope (LF6). Each of these landform classes is identified based on the threshold value of the combination of degrees of slope and curvature described in Table 1 This landform is closely related to the pattern of movement and distribution of water flow, which can explain the morphological properties of the soil of each class of landform elements. In addition, it can also identify the morphoprocess of each landform feature, being erosional on the surface and so can develop and trigger landslides [50].

Land Cover Classification
Geographical Object-Based Image Analysis (GEOBIA) is implemented as the first stage for land cover classification from orthomosaic imagery [51]. GEOBIA is an image segmentation method that groups a set of segments in an image based on a group of pixels that display homogeneous features, namely spectral, radiometric, geometric, and others [51][52][53]. This method is more suitable for very high-resolution imagery, namely orthomosaic from drones, because it can show the presence of massive shadows, low spectral information, and a low signal-to-noise ratio [54,55].
One of the algorithms in the GEOBIA approach is multi-resolution segmentation (Baatz, 2000) that is implemented in the eCognition software to create a set of objects in the image. For segmentation settings, the weights of all three bands (RGB) in the orthomosaic image are equated to be segmented in the scale parameter 150, and the shape/compactness homogeneity criterion is set to 0.3/0.5. In this study, several additional features such as spectral, texture and shape were analyzed through eCognition in each object, followed by feature selection using the Correlation-based Feature Selection (CFS) method in WEKA software [56], as was done in previous research by Ma [56]. Later, the Random Forest (RF) classifier was implemented for land cover classification [57,58] because it is less sensitive to data dimensionality; however, the training sample size was small [57]. In addition, RF is often used as a guided classification for GEOBIA because it can produce land cover and land use maps with good accuracy, both images at medium and very high scales [59][60][61]. Finally, the randomized training and validation samples were used for the RF classifier. This study used sample proportions of about 70% and 30% for training and validation samples, respectively. Classification validation uses a confusion matrix followed by accuracy metrics, namely overall accuracy and kappa coefficient [62,63].

Vegetation Density
The drone only carries a digital camera sensor, which can only photograph objects in the visible spectrum, so it cannot calculate the vegetation density index that requires the near-infrared band. Therefore, vegetation density was analyzed based on the number of vegetation stands per microtopographic unit in km 2 . This is because the vegetation pattern is related to geomorphic processes-including the morphology of the scars-in topographical units [24,25].

Transformation of Categorical Data into Numeric on Microtopography, Landform and Land Cover
Microtopographic zoning, landform, land cover and vegetation density have been constructed to be used as independent variables in the GWR model. However, microtopography, landform and land cover variables are categorical data, which is a problem for the regression model due to data redundancy. To fulfill the requirements in the GWR analysis it is necessary to transform categorical data into numeric. Two data transformation analyses were implemented in this study, namely the proximity factor to feature boundaries based on Euclidean Distance and landscape metrics using FRAGSTATS developed by Dr. Kevin McGarigal with Eduard Ene and Chris Holme as programmers [64]. Thus, this study offers three types of data transformations to find the most suitable GWR model.

Type I Data Transformation
Type I data transformation is based on the distance to feature boundary metrics (using Euclidean Distance) for each class of microtopographic, landform and land cover variables. Previous studies in other fields also used data transformation based on distance metrics, namely land cover variables for spatial modeling [65][66][67]. Technically, each class is analyzed by distance metrics, where the closer to the class feature boundaries, the smaller the distance value (in meters). In addition, the area within the class feature boundary is set to a minus value so that the closer to the midpoint within the class feature area, the greater the distance value with a minus value. This setting is used to distinguish between classes and non-classes and is also considered to represent conventional transformations, namely class 0 (non-class) and 1 (class) dummy variables.
After that, the distance metrics for micro-topography, landform and land cover variables were calculated. Then, a zonal statistical analysis was carried out to calculate the average distance metric for each class that was then aggregated into micro-topographic units.

Type II Data Transformation
Landscape metric analysis is used to transform Type II data, especially the landform and land cover variables as independent variables. This is because the unit of analysis used for the GWR model is the microtopographical unit. In other words, microtopographic units can only be transformed into Type I data. Landscape metrics reflect spatial pattern characteristics. Generally, landscape metrics are often used as predictor variables for ecological analyses, especially for evaluating changes in land cover and land use [68]. However, in this study, landscape metrics-like the distance metric method-were used to transform data into land cover and landform variable categories that could reduce redundancy [69]. In this study, 11 class-level based landscape metrics were taken from several previous studies [69][70][71][72] and implemented to each of the landform and land cover classes, respectively, as shown in Table 2. Landscape metrics analysis utilizes the 'landscape metrics' package (reimplementation of FRAGSTATS) rather than FRAGSTATS software via R language [73]. This is because it can calculate landscape metrics locally and simultaneously (Nowosad, 2022), i.e., per microtopographic unit.

Metrik Rumus Range
Aggregation index (AI) The 11 landscape metrics calculated for each land cover and landform class can produce many features, causing multicollinearity and redundancy between metrics [74,75]. Therefore, the Principal Component Analysis (PCA) approach was implemented to reduce dimensionality by compressing landscape metric features in each class, which was also studied [75][76][77][78][79][80] From the results of PCA calculations, the highest percentage of eigenvalues in the Principal Component (PC) is chosen to represent all landscape metrics for each land cover and landform class, which is named the comprehensive index of landscape metrics.

Type III Data Transformation
The comprehensive index of landscape metrics is also used in transforming Type III data but does not use (PC) with the highest percentage of eigenvalues. Instead, the  2021) developed a combination of PCA and Information Gain methods to reduce dimensionality while selecting the best features of a PC. The result is that it can significantly improve the performance of machine learning classifiers compared to other feature selection methods, including Correlation-based feature selection (CFS), Gain Ratio, and Relief-F. Thus, this study uses the conceptual approach [76] but the aim is to investigate the best PC for each landform and land cover class. In addition, the Relief-F algorithm [77,78] is implemented because it can analyze target data (i.e., landslide data) in numerical form. Relief-F calculates the average merits (AM) in each land cover class and landform which show the ranking of PC attributes. The initial hypothesis for the transformation of Type III data is that it is not always that the Principal Component with the largest percentage of eigenvalues shows a strong contribution to the landslides in the study area.

Geographically Weighted Regression Model Analysis
The GWR model was run to explore the local spatial relationships of landslides with microtopographical variables, landform, land cover, and vegetation density. GWR is a local regression developed by Brunsdon et al. (1998), an update of the Ordinary Least Squares (OLS) method. The GWR model was built based on the percentage of landslide area as the dependent variable and microtopography, landform, land cover, and vegetation density-with all kinds of data transformations carried out-as independent variables that are integrated with microtopographic units. GWR analysis was performed through the GWR4 software originally developed by [79]. The GWR model formula is described as follows [36]: where (µ i , v i ) represents the coordinates of the observed data; β 0 (µ i , v i ) is the intercept parameter at location i; p is the number of independent variables; β k (µ i , v i ) is the local regression coefficient for the independent variable kth at location i; x ik is kth independent variable in ith unit; and ε i is a random error. The regression coefficient is calculated using the local weighted least squares function with the following formula: where W(µ i , v i ) is the spatial weighting matrix of the observation data at the sample point i which represents the effect of sample point i around the regression point on other regression points. In other words, the closer to the sample point i, the greater the influence of the local regression parameter with a larger weight value, and vice versa. Then, the selection of the kernel function is important to determine the scope of the degree of spatial autocorrelation [80]. In this study, the adaptive bi-square kernel was used because the distribution of observational data was not uniform [79]. Here is the bi-square kernel formula: where d ij is the Euclidean distance between sample point j and point i; and θ i(k) is a measure of adaptive bandwidth that shows the spatial variation in the GWR model. Bandwidth selection is also crucial because it measures how well the GWR model generalizes data similar to the data that have been trained. The golden search function determines the optimum bandwidth for the adaptive bi-square kernel function. The optimum bandwidth is determined by the corrected Akaike Information Criterion (AICc) method [79]. The AICc method is known to overcome the problem of over-fitting the model rather than the cross-validation method [80]. When AIC is minimum, bandwidth size is the best.

Results of Mapping Microtopographic Zoning, Landforms, Land Cover, and Vegetation Density
Morphological mapping has been carried out systematically through a remote sensing approach in the study area with a slight modification from the morphological mapping system by Cooke and Doornkamp, namely the Ridge and Valley classes divided into major and minor. Thus, this study's original morphological class numbered 8 was updated to 15 classes (Figure 4). Then, thiessen polygons were applied to construct the microtopographic zoning of each morphological class shown in (Figure 5). A total of 300 microtopographic units were formed in the study area. The symbolization system on micro-topographic maps is based on a combination of colors and textures. Red, purple, blue, brown, and green represent head scarp/cliff, ridge, valley, break of slope, and smooth slope change, respectively. Microtopography is part of the geomorphological study that has an essential role i landslides. Microtopography is defined as changes in topography that can be identifie and mapped on a detailed scale. Based on Figure 4, the microtopography in the study are is divided into four essential parts: (i) ridge, which is divided into the major ridge an minor ridge in the form of sharp and round; (ii) valley, which is divided into major valle and minor valet in sharp and round form; (iii) break of slope which is divided into angula convex, angular concave, and convex and concave too close; and (iv) smooth change i slope, which is divided into convex, concave and convex and concave too close togethe Microtopography is part of the geomorphological study that has an essential role in landslides. Microtopography is defined as changes in topography that can be identified and mapped on a detailed scale. Based on Figure 4, the microtopography in the study area is divided into four essential parts: (i) ridge, which is divided into the major ridge and minor ridge in the form of sharp and round; (ii) valley, which is divided into major valley and minor valet in sharp and round form; (iii) break of slope which is divided into angular convex, angular concave, and convex and concave too close; and (iv) smooth change in slope, which is divided into convex, concave and convex and concave too close together. Each form of microtopography has a different effect on landslides. The movement of soil material will increase on sharp slope morphology. In addition, differences in microtopography will also affect the value of the shape of the land surface in the form of plan curvature and curvature profile. Both influence the acceleration and deceleration of the water flow. In ridge microtopography, the water flow will be accelerated so that the potential for material movement increases. The pattern of microtopographic distribution in the study area can be said to be heterogeneous. The southwestern part is dominated by ridge and valley morphology with various class variations. The shape of the zoning in this area tends to be elongated, reflecting the presence of hills. The major valley (sharp) with a large proportion of the area is a perennial river channel stretching east to west. Then, the central part is formed by a fragmented microtopography, namely the form of zoning with a relatively small proportion of areas and various morphological classes. This implies complex slope configurations and slight fluctuations. The northeastern part is almost the same as the southwest, namely the form of an elongated zonation, but the morphology is relatively fragmented. Mean- The pattern of microtopographic distribution in the study area can be said to be heterogeneous. The southwestern part is dominated by ridge and valley morphology with various class variations. The shape of the zoning in this area tends to be elongated, reflecting the presence of hills. The major valley (sharp) with a large proportion of the area is a perennial river channel stretching east to west. Then, the central part is formed by a fragmented microtopography, namely the form of zoning with a relatively small proportion of areas and various morphological classes. This implies complex slope configurations and slight fluctuations. The northeastern part is almost the same as the southwest, namely the form of an elongated zonation, but the morphology is relatively fragmented. Meanwhile, the northern part shows a more compact, elongated form of microtopographic zoning and a large proportion of areas indicating a less volatile slope configuration.
Then, the shape-based landform variables by the Pennock system have been estimated and classified. In contrast to microtopographic zoning, the Pennock landform in the study area is defined by curvature and slope indices, namely profiles and plans with a neighborhood radius of 40 m each. As shown in Figure 4a, the predetermined neighborhood radius scale reflects the detailed landform pattern. Generally, in the Pennock method of landform the slope is divided into three arrangements, covering the upper (shoulder), middle (backslope), and lower (footslope), with almost the same proportions. However, the Pennock method classifies a convergent/divergent backslope of around 3% of all existing landforms in the study area. Meanwhile, the proportion of convergent/divergent shoulders and footslopes is almost equal, but the spatial pattern still varies. In addition, the landform in the study area tends to have a convex-concave pattern that repeats over short distances.
Class 1 hierarchical land cover classification on orthomosaic imagery based on GEOBIA shows very good accuracy, with an overall accuracy of 98.26% and a kappa coefficient of 0.97. Table 3 presents the results of the confusion matrix calculation using the Random Forest guided classification. The selection of attributes on spectral features, texture, and shape index also influences satisfactory accuracy. Of the 27 features used as predictor variables for classification, 13 features were selected based on the results of attribute filtering using CFS (Correlation-based Feature) from WEKA software. In Figure 4b, most of theland cover in the study area is dominated by herbaceous (45.77%) and forest (36.56%), followed by bare soil (12.71%), shrub (3.58%), and built-up (1.36%). In addition, the density of moderate to very high vegetation classes has spread from the north to northeast and slightly in the center and southwest (Figure 6c).

Comprehensive Index Analysis of Landscape Metrics on Landform and Land Cover
A comprehensive index of landscape metrics for each land cover and landform class is analyzed using the PCA algorithm to compress metric landscape features. As shown in Tables 4 and 5, the highest percentage of eigenvalues is PC1, which is overall above 30% and the eigenvalues are above 1. This indicates that the PC1 component has most of the feature information from all landscape metrics, so it is used for data GWR Type II model independent variable input. In addition, while each PC is associated with landslide data and the results are not always PC1, it has a large contribution to explaining landslides in the study area, though best represents all landscape metrics. Interestingly, the LF1 and LF6 classes show that PC6 contributes the most, although the eigenvalue is below 1. In fact, negative AM values appear in PC1, namely classes LC2 and LC5, which means that it has the lowest contribution among other classes to landslides, even though the negative AM value is still used for the GWR model input.
dictor variables for classification, 13 features were selected based on the results of attrib filtering using CFS (Correlation-based Feature) from WEKA software. In Figure 4b, m of theland cover in the study area is dominated by herbaceous (45.77%) and fo (36.56%), followed by bare soil (12.71%), shrub (3.58%), and built-up (1.36%). In addit the density of moderate to very high vegetation classes has spread from the north to no east and slightly in the center and southwest (Figure 6c).

Model the Local Spatial Relationship between the Landslide and the Selected Independent Variable Model Type for Priority Land Management Sustainability
In this study, the GWR model was implemented to investigate the spatial pattern of local landslide relationships with independent variables collected in spatial units of microtopographic areas. Table 6 shows the GWR fit model on the data group of predictor variables influencing landslide events (microtopography, landform, land cover, and vegetation stand density), adjusted for R 2 and AICc. All models use the Gaussian kernel, whose bandwidth size is checked and calculated based on the Golden Search method. AICc is used as a selection criterion to find the optimal bandwidth of 300 observational data.
Overall, in the univariate local sub-model the LC model Type I variable group indicates the covariate most related to the landslide event, namely the value of adj. The highest R 2 (0.789) and the lowest AICc (2478.33) are even more important than all the univariate and multivariate sub-models analyzed. A multivariate model of two groups and three groups of variables, the Type I model shows better performance (adj. R 2 above 0.7) compared to the Type II and Type II models. The univariate and multivariate sub-models, based on the type of model, suggest that Type I shows the best fit model compared to Type II and Type III. In addition, the Type III model is better than the Type II but performs lower than Type I. On the other hand, the sub-model with all variable groups included shows a different pattern of Type model results, namely Type III is the highest (adj R 2 = 0.755) with the third lowest AICc value (2526.38), after the LC group Type III (2515.51) and Type I (2478.33) sub-models. Thus, the Type III group sub-model with four variables was chosen to provide additional insight, namely the local spatial relationship experiment of the landslide. However, in terms of calibration the fit model is not the best compared to other sub-models. Table 6 is a summary of the results of univariate and multivariate comparisons of the GWR model with four groups of predictor variables: microtopography (MT); landform (LF); land cover (LC); and standing vegetation density (VD). Each predictor variable has subvariables for each class whose data form is divided into three types: proximity factor (Type I); PC with the highest percentage of eigenvalues from the landscape metric comprehensive index (Type II); and PC on the comprehensive metric index.
In the summary output of the GWR model (Table 6), the F-statistics in the ANOVA comparison test shows that the entire GWR model significantly increases global model performance (OLS). Thus, the null hypothesis of the GWR model being unable to improve the performance of the global model is rejected. Figure 5 shows the local R 2 generated from the GWR, that the GWR model is fit to map local landslides, which are explained through the MT, LF, LC, and VD predictor variable groups that are modeled simultaneously. The distribution of R 2 in the study area tends to be homogeneous and clustered. About 80% of the local model area has an R 2 above 0.73, and at least 10% of the area in the southwestern region shows that the local model is less fit (fit) with an R 2 below 0.43. This implies that additional covariates are needed to explain the slides in the study area, particularly the southwest region. Then, the positive value of the standardized residual indicates overestimation. Likewise, negative values that appear indicate an underestimation. Overall, the GWR model shows an over-/under-estimated distribution pattern that tends to be random in the study area (Figure 7). The local coefficients shown in Table 7 indicate that the relationship between the landslide and the independent variables is non-stationary. The relationship varies spatially with a range of magnitudes and directions. The independent variables in the MT, LF, LC, and VD groups showed that the magnitude of the relationship varied in both negative and positive directions in terms of min and max values. This can be interpreted as the presence of an increasing variable that will also increase the occurrence of landslides, but the presence of variables can also reduce the occurrence of landslides. However, several independent variables only show an inverse correlation to landslides, including Micro2, Micro6, LF3, and LC3. This suggests that these four variables do not contribute to landslides within this feature area. In addition, the variables LF4 and LF6 positively contribute to the occurrence of landslides in all study areas in this feature. However, the VD variable has the smallest relationship magnitude compared to the other variables, with a two-way relationship.
comparison test shows that the entire GWR model significantly increases global model performance (OLS). Thus, the null hypothesis of the GWR model being unable to improve the performance of the global model is rejected. Figure 5 shows the local R2 generated from the GWR, that the GWR model is fit to map local landslides, which are explained through the MT, LF, LC, and VD predictor variable groups that are modeled simultaneously. The distribution of R2 in the study area tends to be homogeneous and clustered. About 80% of the local model area has an R2 above 0.73, and at least 10% of the area in the southwestern region shows that the local model is less fit (fit) with an R2 below 0.43. This implies that additional covariates are needed to explain the slides in the study area, particularly the southwest region. Then, the positive value of the standardized residual indicates over-estimation. Likewise, negative values that appear indicate an underestimation.
Overall, the GWR model shows an over-/under-estimated distribution pattern that tends to be random in the study area (Figure 7). The local coefficients shown in Table 7 indicate that the relationship between the landslide and the independent variables is non-stationary. The relationship varies spatially with a range of magnitudes and directions. The independent variables in the MT, LF, LC, and VD groups showed that the magnitude of the relationship varied in both   Generally, in local relationship analysis research the local coefficients of each independent variable in the GWR model are imported and visualized as a map using the GIS environment. The coefficient values are mapped taking into account the two-tailed t value, i.e., t values above 1.96 and −1.96 are considered significant (equal to p < 0.05). However, in this study the local coefficient maps are summarized in one map by integrating the four groups of variables based on the largest local coefficient values (negative or positive values). This is because the spatial units used are microtopographic units, so the other three variables are superimposed and the value of the feature slices included in each unit is calculated. In this case, the value of the principal component comprehensive index of the selected landscape metrics is calculated. The estimated values of parameters or local coefficients on the microtopography are searched for and adjusted to each class of unit so that misinterpretation does not occur, e.g., extracting the Micro1 coefficient values based on the attributes of the Micro1 class microtopographic units so that outside the area has been selected. Likewise, with the landform and land cover, only the unit containing information from the two variables is taken for the coefficient value.
As shown in Figure 8, there are 14 classes combined with the local coefficients of the independent variables of each group. Signs "+" and "−" are interpreted as the relationship's direction. Map visualization uses a combination of colors and textures to make it easier to read map symbols and their information. Each unit of the analysis found several variables positively or negatively related to landslides. In other words, the characteristics and patterns of landslides in the study area can be explored through several variables related to the unit of analysis.
In the map in Figure 8, MT classes that show a positive relationship to landslides are green, and tend to be spread over the central area, being most prominent in the northeast. Class LF2, LF6, and LC2 also show a positive relationship. However, units with no significant positive or negative relationship are clustered in the southwestern region with white symbols. This can also be attributed to the low R 2 value in the region, as shown in Figure 7. Interestingly, the units that only show a negative relationship, namely the MT, LC3, and LC4 classes, appear in the northern region with a grayish color symbol. Vice versa, which only shows a positive relationship appears in the south with a dark purple color, but only one unit, namely the LF6 class.
As shown in Table 8 Meanwhile, only micro2 (1 unit) and micro6 (3 units) have a negative relationship. Thus, LF6 and LC3 have an important role in understanding the landslide mechanism in the study area. Figure 9 shows a map inset focused on landslide and non-slip areas overlapping the most significant micro-topographical units and classes. The current landslides were associated with LF6, rather than MT and LC2. Regardless of the significance of the local relationship, the current landslides predominately occur at micro6, micro8, micro9, and micro13.
the attributes of the Micro1 class microtopographic units so that outside the area has been selected. Likewise, with the landform and land cover, only the unit containing information from the two variables is taken for the coefficient value.
As shown in Figure 8, there are 14 classes combined with the local coefficients of the independent variables of each group. Signs "+" and "−" are interpreted as the relationship's direction. Map visualization uses a combination of colors and textures to make it easier to read map symbols and their information. Each unit of the analysis found several variables positively or negatively related to landslides. In other words, the characteristics and patterns of landslides in the study area can be explored through several variables related to the unit of analysis. In the map in Figure 8 , MT classes that show a positive relationship to landslides are green, and tend to be spread over the central area, being most prominent in the northeast. Class LF2, LF6, and LC2 also show a positive relationship. However, units with no significant positive or negative relationship are clustered in the southwestern region with white symbols. This can also be attributed to the low R 2 value in the region, as shown in Figure  7. Interestingly, the units that only show a negative relationship, namely the MT, LC3, and LC4 classes, appear in the northern region with a grayish color symbol. Vice versa, which only shows a positive relationship appears in the south with a dark purple color, but only one unit, namely the LF6 class.

Mapping Microtopography, Landforms, Land Cover, and Vegetation Density
Microtopographic zoning mapping can be said to be rarely discussed in the field of geomorphology. In addition, Yang et al. used the Thiessen polygon approach to construct karst landform zoning precisely based on the spatial proximity of peak and nadir points, namely positive and negative landscapes [81]. In this study, microtopographic zoning was also derived from the results of thiessen polygon-based buffering on the morphology of the Cooke and Doornkamp system. Each mapped morphological boundary is used as a Thiessen polygon control point to build micro-topographic zoning. The results show that

Mapping Microtopography, Landforms, Land Cover, and Vegetation Density
Microtopographic zoning mapping can be said to be rarely discussed in the field of geomorphology. In addition, Yang et al. used the Thiessen polygon approach to construct karst landform zoning precisely based on the spatial proximity of peak and nadir points, namely positive and negative landscapes [81]. In this study, microtopographic zoning was also derived from the results of thiessen polygon-based buffering on the morphology of the Cooke and Doornkamp system. Each mapped morphological boundary is used as a Thiessen polygon control point to build micro-topographic zoning. The results show that the shape and size of the microtopographic units vary ( Figure 5). The large and rounded unit shapes indicate that the distance between morphological boundaries is farther apart than the small and elongated unit shapes. It can also reflect the degree of geomorphic processes to form different microtopography. In other words, different morphologies allow different morpho processes to form different surface material characteristics [10].
Elemental landforms based on the Pennock system have been classified using highresolution DTM data. The radius scale for the plan and profile curvature, which is 40 m each, is set based on more detailed spatial variability, but this is visual, so there may be uncertainty in the classification. In addition, the Pennock method identified the slope arrangement of the shoulder and footslope sections on both convergent and divergent slope shapes. Unfortunately, the Pennock method seems less sensitive for identifying the backslope. This problem is contrary to the study of Evans et al. that the Pennock method failed to identify the arrangement of shoulder and footslope slopes on agricultural land with low relief configurations [49]; it is also reported that the study area could not map divergent footslope landforms but was dominated by divergent back slopes with more complex relief configurations [82]. Thus, the differences in the problems in this study with other studies may be caused by differences in spatial resolution, unique relief configurations, and setting the radius parameter scale for which there is no definitive rule.
GEOBIA works well in class 1 hierarchical land cover classification with orthomosaic imagery: forest, shrubs, herbaceous plants, open land, and built-up land. Overall, producer accuracy and user accuracy for all classes reach above 90%. For the user, accuracy metric value of built land is the lowest among the other classes. This may be due to several built-up land objects having similarities in spectra, shape, and texture with forest areas, herbaceous plants and open land. However, this problem is not significant because it basically meets the land cover mapping requirements, which are above 85% [83].
There are several reasons why the accuracy of this land cover map is satisfactory. First, the method offered by [60], namely the combination of GEOBIA and selection of spectral features, shapes and textures, can classify land cover from orthomosaic imagery very well. Furthermore, segmentation parameters such as scale, shape, and compactness are obtained through trial and error, especially for scales set with a range of 50 to 500 in multiples of 50. As explained in the method Section 2.3.4., the optimum scale parameter is found, namely 150. Third, the CFS method works well in selecting all features into 13 optimal features for land cover classification in the study area. In addition, the RF classifier also performs very well, which is even slightly superior to the study by De Luca et al. with orthomosaic imagery equipped with near-infrared bands [61].

GWR Model Implementation with Different Independent Variable Data Transformation
Categorical data transformation for each class of independent variable has been carried out for the needs of the GWR model. The Type I model represents numerical data based on distance metrics microtopographic, landform and land cover variables. Then, the Type II and Type III models represent a comprehensive index of landscape metrics for each class of land cover and landform variables using PCA. These three types of models are crucial for the GWR model to work well because of the conventional-based data transformation problems, namely binary classes. However, there is a uniqueness in the Type III model, which breaks the standardization of PCA analysis that PC1 retains most of the information by maximizing the variance of the data from the comprehensive index of landscape metrics compiled in the Type II model. In other words, the PC with the highest AM value has the greatest correlation with landslides, although the eigenvalue is less than 1. Class variables with low variation are not unimportant [84] when associated with landslide data. This case is like the climate study by Jolliffe, in which the low variation component relates to the response variable rather than the high variation component [85]. In addition, discarding PCs with small eigenvalues can lead to bias [86]. Hadi and Ling also reported that PCs selected based on the breakdown of principal components that depend only on the variation in variable X (i.e., the independent variable) might fail to account for regression fit because they do not contribute anything to the response variable (or the dependent variable) [87]. It was proven in the results of this study that when the Type II and Type III models were analyzed byboth univariate and multivariate GWR, the overall values of AICc and adj. R 2 in Type III is higher than Type II (Table 6). In addition, when the Type III model on landform and land cover is combined with the microtopography and vegetation density variables that are modeled simultaneously, it shows a higher regression fit than the variables in the Type I and Type II models. Meanwhile, the univariate model in Type I showed a higher regression fit than Type II and Type III. This indicates that setting compression and selecting landscape metric features on landform and land cover can increase the predictive power of the GWR model for exploring local spatial relationships of landslides in the study area.

Modeling Local Spatial Relationships of Landslides
The GWR model has been analyzed to identify the local spatial relationship between the landslide and the micro-topography, landform, land cover, and vegetation density in the study area. The Type III model with the best regression fit among the others was chosen as the independent variable input and the percentage of landslide area as the response variable. The result is that most of the independent variables can explain the landslides in the study area but not for the southwestern region. The statistical summary of parameter estimation (Table 7) shows that almost all variables spatially have positive and negative estimates which indicate a non-stationary relationship to landslides. The local coefficient integration map is the result of a synthesis of all variables that have been selected based on their significance and adjustment of feature information for each variable (namely, local coefficients that intersect with their features). This map is more intuitive and less confusing in examining the most significant variables, both the positive and negative relationship direction to landslides in each unit ( Figure 6).
Variables positively related to landslides include microtopography, convergent backslope, divergent footslope and shrubs. In contrast, those negatively related include a small portion of microtopography, herbaceous plants, and open land. Meanwhile, the variables that showed the most positive and negative relationships were divergent footslopes and herbaceous plants, respectively. Sato et al. reported that landslides in the mountains of the northwestern Himalayas occurred relatively on convex slopes rather than concave slopes [88]. The study by [89] reported that landslides in Changshou Valley, Baoji City in Shaanxi Province often occur in surface relief with convex slope shapes. Pourghasemi et al. also reported that the shape of the convex slope had a major effect on landslides in the Jumunjin Area, South Korea [28]. Havenith et al., revealed that the reasons why landslides occur on convex slopes include: (1) convex slopes are relatively less stable under similar hydrogeological conditions (lower factor of safety) because a larger slope body (larger driving force) acts on the same sliding surface (more the same resistance); and (2) they allow for the presence of a accumulated material (colluvium), which reflects lower shear strength [90].
This statement supports the results of this study, namely that convex slopes are closely related to landslides in the study area. In addition, this convex slope is specifically located at the foot of the slope with surface material that holds a higher water content than the backslope and upper slope arrangement [48]. Then, when the water content in the surface material increases, the pore water pressure increases, which can result in a low Factor of Safety (increased shear stress and decreased shear strength) resulting in landslides [91]. In addition, the dynamics of surface material moisture content can be affected by rainfall intensity. As additional information, the divergent footslope in the study area is under an average slope of 28.8 • ± 15.2 • , thereby increasing the potential for landslides to occur. Specifically, herbaceous plant land cover is plantation land that is negatively related to landslides in the study area. This finding can be explained as land that was not initially maintained, such as shrubs and grasslands; open land converted to cultivation implies better water and land management practices that can allow for reduced slope instability [27].
However, several other significant variables also need not be ignored; although, only a few units.
Based on the findings of the spatial relationship of landslides with micro-topographical and vegetation variables, it can be used as a basis for the conservation of research areas to reduce the potential for landslides. Based on the results of the GWR, it shows that landslides occur on convergent backslopes with surface cover in the form of shrubs. In addition, landslides also have the potential to occur on convex slopes that are on the toe slopes. These findings can become a model for micro-topography and vegetation-based conservation arrangements, so that to reduce landslides the community in the study area is not advised to plant species in the form of shrubs on steep slopes, or to carry out intensive processing on footslopes with convex slope shapes. Herbaceous plants show a negative relationship to landslides, so as part of the landslide conservation effort planting of herbaceous plant species such as Rumput gajah (Pennisettrum Purpureum), Pakis (Diplazium esculentum), Mindih (Melia azedarach L), and Waru (Hibiscus tiliaceus) can be applied. Several herbaceous plants in the study area are conservative if planted using multi strata techniques, thereby reducing the potential for movement of soil material.

Research Implications and Limitations
This study reveals local relationships between landslides and microtopography, landforms, land cover, and vegetation density. However, the novelty offered is to build a natural unit of analysis compared to a grid basis, which is based on micro-topography derived from morphological thiessen polygons. Pennock system-based landform element independent variables are also taken into account. To our knowledge, they are rarely analyzed for detailed scale landslide studies because they are indirectly related to geomorphic processes. Land cover classification uses a reliable method that produces the best accuracy-using the method from the study of Ma et al. applied to drone-based orthomosaic imagery, which can be replicated for other studies. Then, three types of transformation of categorical data on independent variables become new insights with an analytical approach for the study of regression-based spatial modeling that can be reapplied, especially in research on geo-environmental disasters. In addition, the last and most important thing is to offer the output of the GWR model with a local coefficient integration map for the most significant classification of variables related to landslides (positive or negative). This map information can be used as a reference for landslide disaster management in the study area, namely convergent backslope, divergent footslope, some microtopographical classes, and shrubs that need to be watched out for in their land use and require special treatment to reduce slope instability.
However, there are several limitations in this study, including: (1) the mapped study area contains canopied vegetation, which is a weakness in building DTMs, so that it can cause land surface elevation estimation errors; (2) microtopographic zoning does not have field validation, so objectivity is still questionable because this study relies on interpretation remote sensing; and (3) the results of the GWR model with the transformation of landscape metric-based category data are not discussed further as to how the spatial configuration of landforms and land cover can specifically explain landslides in the study area. For further research, additional issues such as climate and anthropogenic factors can be added to the GWR model to obtain more complete results and to explain landslides in the entire study area.

Conclusions
Exploration of local spatial relationships between landslide occurrences in tropical hills, especially in Taji Village, Jabung District, East Java Province, can explain the factors that most influence landslides. A series of approaches and methods were introduced and implemented to construct independent variables to be analyzed by GWR and produce reliable model outputs. Microtopographical zoning is used as a unit of analysis to synthesize all information on the independent variables-which incidentally have different spatial forms-as well as the dependent variable, namely the percentage of avalanche area. Categorical data transformation brought independent variables to be modeled, even better than global regression. However, the best results for the independent variables modeled simultaneously fell to the change in landscape metric data on landform and land cover features that had been compressed through PCA analysis, and selected component features associated with landslide occurrence data. GWR can reveal a non-stationary relationship between landslide events and independent variables in the study area. Information on the local coefficients of each independent variable is integrated into a single entity based on the significance of the t values (i.e., ≤−1.96 and ≥1.96) and the selection of features of each variable that intersect with the microtopographic unit. The majority of the variables that show the most positive relationship to the occurrence of landslides are divergent footslopes where colluvial (colluvial) material accumulates from the upper slopes (i.e., the shoulders and backslope), which have a low Factor of Safety so that the slopes are less stable. In addition, the slopes of the footslope accommodate more water content, and when heavy rain occurs it can increase the pore water pressure so that the Factor of Safety is lower and there is a potential for landslides to happen, especially since the angle of the slopes in the study area is sufficient to support this. On the other hand, it was herbaceous plants or plantation land in the study area that surprisingly reduced the occurrence of landslides, due to good water and land management being able to maintain slope stability longer. The methodological approach developed and introduced in this study is reproducible and further analyzed in other tropical hills at a detailed and regional scale prone to landslides.