The study was conducted within the ecological boundary of the Opak Oyo Watershed (OOW) (Figure 1), which is administratively located in Bantul Regency, Sleman Regency, Gunungkidul Regency and Yogyakarta City, covering an area of approximately 1304.95 km². The boundaries of the catchment area were identified using data from the Ministry of Environment and Forestry. Topographically, the area extends from the upper volcanic slopes of Mount Merapi in the north to the karst-dominated landscapes of Gunungkidul in the south, with elevations ranging from approximately 50 m to over 2,900 m above sea level. Slope gradients vary considerably, from flat alluvial plains (<8%) in the downstream areas, to steep volcanic terrains (>40%) in the upstream zone.

Land cover within the watershed is heterogeneous, consisting of protected forest in the upper slopes, mixed agriculture and plantations in midstream areas, and expanding built-up zones in peri-urban and tourism corridors. Forested areas play a critical role in slope stabilization and sediment retention, while agricultural and tourism-driven land conversion has intensified spatial pressure in recent decades. These physical and environmental characteristics make the OOW a highly dynamic volcanic landscape, where geomorphological processes, hydrological dynamics, and tourism development interact simultaneously.

The research was conducted over three main stages: pre-field preparation, field verification, and post-field data processing. The pre-field stage focused on collecting and preparing secondary spatial datasets and analytical variables. The compiled variables included multi-temporal Landsat imagery (2004, 2014, 2024) for land-use and land-cover (LULC) analysis, DEMNAS-derived topographic parameters (elevation and slope gradient), lithology maps, river networks, road networks, administrative boundaries, population density, and official Volcano Hazard Zone (KRB) maps. Driving factors used in the Land Change Modeler (LCM) consisted of distance to rivers, roads, economic activity centers and government centers.

The field stage was conducted to validate the land-use classification results and tourism site locations through ground truth surveys using GPS coordinate recording and direct observation of land-cover categories (forest, agriculture, built-up areas, water bodies and open land). Tourism facilities were verified in terms of their spatial position, surrounding land use, and proximity to hazard zones. The post-field stage involved classification accuracy assessment using confusion matrix analysis, LCM-based land-use change simulation up to the year 2034, and spatial overlay analysis between tourism sites, topographic layers and KRB maps. The final outputs included thematic maps, statistical summaries and spatial interpretation to support sustainable tourism planning. The overall workflow distinguishing between pre-field, field and post-field activities is presented in Figure 2.

Secondary data, such as road data, river data, lithology data and DEMNAS data, were used, obtained from several institutions. The secondary data were then processed utilizing remote sensing and GIS integration to create topographic maps and land-use change prediction maps in the OOW. The data employed can be seen in Table 1.

The Earth's physical surface conditions can be represented by topographic maps showing elevation and slope. Elevation refers to the height of a location above sea level, while slope refers to the percentage ratio between vertical distance (land height) and horizontal distance (flat land length) ((Darmawan et al., 2017); (Triwahyuni, 2017)). Slope shapes are generally influenced by erosion, land movement, and weathering processes. Information regarding the topography of the OOW was obtained by utilizing the National Digital Elevation Model (DEMNAS) with geospatial technology, as shown in Figure 4.

The topographic map of the OOW (Figure 4) shows significant variations in elevation and slope gradients. Based on the elevation map (Figure 4a), the area has an elevation range from 0 to over 300 meters above sea level. The southern part of the OOW is dominated by low-elevation areas (0–50 meters), while the northern part has higher elevations, reaching over 300 meters, which are likely to be hilly or mountainous regions. Areas with higher elevations generally have the potential to function as water catchment zones, whereas lowland areas are more vulnerable to surface runoff and flood risks. Table 4 shows the area coverage for the elevation in the Opak Oyo Watershed. Meanwhile, the slope gradient (Figure 4b) indicates that the OOW has a diverse topography, ranging from flat plains to very steep hills. Table 5 shows the area coverage for slope gradients in the watershed.

In the 3D visualization (Figure 5), it was observed that the central part of the OOW is bounded by hills with steep slopes. Furthermore, the region is dominated by lowland areas with gentle slopes, gradually descending towards the downstream areas, which directly border the sea.

The CART algorithm is widely used to analyze land-use changes ((Epie & Hull, 2025); (Tekla et al., XXXX)). Based on the data processing conducted using Google Earth Engine, it was found that the results were reliable and could be used for further processing, as shown in Table 6.

A sampling method involving the stratified random sampling technique was used to determine the location points of the previously processed 2024 land use map results. This was done to validate the results of the satellite image classification processing. This method was employed so that the sample distribution was evenly distributed across all land use classes (all classes were represented). The maps processed using the CART algorithm with the training data used consisted of 582 points for total landuse, which were then validated using Google Earth imagery with 182 validation points, resulting in a user accuracy of 0.945, a kappa accuracy of 0.928, and an overall accuracy value of 0.945 (or 94.5%).This indicated that the map was accurate and could be considered valid for further use, as its accuracy was >0.80. Based on the results of the data processing, a temporal analysis was then conducted on the land-use changes in the OOW area. It was found that from 2004 to 2024 significant changes had taken place (Table 7).

Based on Table 7, during the period from 2014 to 2024 there was a significant increase in built-up areas, with the addition of 47.12 km². This was part of the local government's efforts in rehabilitation, reconstruction and recovery in the housing, infrastructure, social and productive economic sectors as a result of the eruption of Mount Merapi and the earthquakes in the region. Furthermore, land use in the OOW was shown to be generally characterized by land conversion from dryland agriculture, shrubs, forests and open land into built-up areas. This change in land-use patterns led to a reduction in water absorption capacity, increasing surface runoff (Nugroho et al., 2018). A visualization of land-use changes in the OOW from 2004 to 2024 can be seen in Figure 6.

Water, agricultural and open land types also saw increases in area. The expansion of water areas from 2004 to 2024 was part of an integrated watershed management effort, particularly for the OOW, as stipulated by the Governor's Decree of the Special Region of Yogyakarta No. 285/KEP/2014. The increase in agricultural areas during the period 2004-2024 in the region influenced the hydrological response of the Opak Watershed. Agricultural land dominated the Opak Watershed, covering 518.36 km² in 2024. Generally, the increase in agricultural land use, particularly for mixed dryland agriculture, was balanced by the expansion of settlement areas. Agricultural land in the OOW was utilized for the cultivation of food crops such as corn, upland rice, soybeans and cassava. According to a study simulating existing land use dominated by agriculture and settlements (Widiatmoko et al., 2020), the KAT (Annual Runoff Coefficient) value was high at 0.4, with a River Regime Coefficient of 85.16, also in the high category. This suggests that a significant amount of rainfall became direct runoff. These findings indicate that land-use changes impacted the hydrological response of the watershed. Therefore, a watershed management planning strategy is needed to optimize the hydrological function of the area. Additionally, an increase in open land area of approximately 5.03 km² occurred between 2004 and 2024.

Land use in the OOW in 2004 and 2014 was employed to develop a Land Change Modeler (LCM) model, with a scenario of increased built-up land in 2024, using the MLP Neural Network-Markov Chain algorithm. In our study, the algorithm was used with the following criteria: a) a final learning rate of 0.0001; b) a momentum factor of 0.5; c) a sigmoid constant of 1.0; and d) eight input layer neurons, four hidden layer neurons and two output layers. These criteria were based on RMS of 0.01, 10000 iterations, and an accuracy rate of 100%. An input layer of eight neurons was used to represent the number of driving variables used in the model, while the hidden layer was employed to process complex non-linear relationships.

The low learning rate value of 0.0001 indicates that the model learns very carefully and slowly, which is good for achieving high accuracy without overshooting the optimal solution, although it requires longer computation time. A momentum factor value of 0.5 is balanced between maintaining the previous learning direction and adjusting to new data. An iteration value of 10,000 indicates the maximum number of learning cycles (epochs), while an RMS error value of 0.01 is the average error tolerance limit. If the average prediction error falls below 0.01, the model is considered sufficiently accurate and stops learning. The LCM simulation with a built-up land increase scenario of OOW in 2024 was built on a standard but rigorous algorithmic foundation. The validity of the prediction results (Table 8) depends heavily on how well the eight input variables (neurons) represent the reality of land change drivers in the field, as the machine parameters themselves are set for high mathematical accuracy. These criteria represent the default settings of the LCM model. The results are shown in Table 8.

The results of the calibration and validation of the LCM model using the multi-layer perceptron (MLP) algorithm show varying but robust levels of reliability across all land transition variables. These findings are in line with recent studies that consider MLP to be one of the most effective algorithms for handling the non-linear dynamics of land cover change compared to conventional methods (Subiyanto et al., 2019). As detailed in Table 8, the model shows high stability in the learning process, characterized by consistently low training RMS and testing RMS values, with a very small difference between them. For example, in the transition from vacant land to built-up land, the training RMS was recorded at 0.3751 and the testing RMS at 0.3744. The closeness of these values confirms that the model experienced no overfitting and was able to generalize well on the test data.

The close RMS values in the query indicate that the model was effectively validated and could generalize well to new data, a crucial indicator of validity in neural network-based spatial modeling (Sabiri et al., 2022). Overall, the highest accuracy rate was achieved in predicting the change from open to built-up land, at 90.44%, with a very dominant transition skill (0.8901). This shows that the MLP algorithm is very effective in capturing spatial patterns in open land. Although the accuracy for the transition to agriculture is slightly lower (67.25%), the model still maintains acceptable predictive capabilities given the complexity and heterogeneity of agricultural patterns, which often pose challenges in tropical region modeling (Sajan et al., 2022). Based on the above MLP results, it was found that the ROC (receiver operating characteristic) value ranged from 67.25% to 90.44%, falling between 50% and 100% (Table 8), indicating that the model used, along with the eight variables, was significant and influenced land-use changes. As demonstrated in Table 8, it was found that each model had dominant and non-dominant factors affecting land-use changes in the Opak Oyo Watershed.

The land-use modeling results using the Land Change Modeler MLP-Markov Chain model with a built-up land expansion scenario for 2024 (Figure 7b), were compared with the actual land-use map for 2024 in the OOW (Figure 7a). Kappa statistics for quantity and location were obtained based on an evaluation of the predicted compared to the observed land use. The statistics show that the Kno value was 0.8176, the Klocation value was 0.8053, the KlocationStrata value was 0.8053, and the Kstandard value was 0.7206. All kappa index values obtained between 0.72 and 0.80 indicate significant agreement, with values between 0.81 and 0.99 indicating near-perfect agreement (Azizah et al., 2023).

Furthermore, the validated model for 2024 was used to model land use with a built-up land expansion scenario for 2034. From the predicted land use map for 2034 (Figure 8), it was found that in that year (Table 9), the area of water bodies in the OOW will cover 22.86 km², forest areas will cover 523.94 km², built-up land will cover 213.87 km², agricultural areas will cover 422.12 km², barren land will cover 5.55 km², open land will cover 4.92 km², and cloud cover will account for 100.15 km². Compared to the 2024 land-use results, several land-use category experienced significant changes and area expansions. Reductions occurred in water bodies, agricultural areas, barren land, open land and cloud cover. In contrast, other categories, such as forests and built-up land, experienced area expansions, at 6.95% for forests and 3.53% for built-up land, between 2024 and 2034. However, water bodies decreased by 0.12%, agricultural areas by 7.67%, barren land by 0.56%, open land by 0.10%, and cloud cover by 2.04% between 2024 and 2034.

Temporal analysis of the changes in land cover revealed two contrasting transition patterns between the historical observation period (2014–2024) and the prediction scenario (2024–2034). During the decade from 2014 to 2024, the landscape of the region was dominated by significant deforestation, with forest cover shrinking dramatically by 127.76 km² (-22.7%). This decline correlates positively with the expansion of agricultural areas, which increased by 59.43 km², and of built-up land, which increased by 47.12 km². This indicates strong anthropogenic pressure, with the conversion of forest land to cultivation and settlement areas being the main drivers of landscape change.

However, the 2034 land use scenario projects a reversal of ecological trends. The agricultural sector is predicted to experience a massive contraction, with a reduction in area of 100.14 km². This decline appears to be offset by land rehabilitation programs or natural succession, marked by the restoration of 90.75 km² of forest areas. Several forest and land rehabilitation initiatives are being undertaken, such as implementation of the Indonesia's FOLU Net Sink 2030 program by the Ministry of Environment and Forestry, and the Tree Collective Gunungkidul 2025 program, which involves the planting of 15,000 multi-purpose tree species (MPTS). This program is organized by the Javlec Indonesia Foundation, JejakIn, and the Wana Makmur and Wana Lestri Community Forest Farmer Groups (Dharmawan & Pratiwi, 2023), and involves the conversion of abandoned or critically degraded Sultan Ground (SG) and Pakualaman Ground (PAG) around the OOW into conservation forests or green open spaces.

Despite fluctuations growth in the forestry and agriculture sectors, urbanization shows a consistent linear growth trend. Built-up land is projected to continue increasing by 46.04 km² during the 2024–2034 period, absorbing vacant land and some agricultural land, bringing the total built-up area to a peak of 228.13 km². This phenomenon confirms that even though ecological restoration efforts are projected to occur, physical urban expansion remains a permanent variable in the spatial dynamics of the region.

In Figure 8, it can be seen that the conversion of land into built-up areas has primarily spread from the central part of the OOW toward the coastal areas. This is shown by the projected potential for transition, where the red visualization with a value of 1 indicates the potential for land conversion into built-up areas in the OOW (Figure 8). As shown in Table 8, the most dominant factor influencing the conversion of water bodies into built-up land is the distance from rivers; the conversion of forests and agricultural areas into built-up land is influenced by the distance from roads; the conversion of barren land into built-up land is influenced by the distance from economic activity centers; and the conversion of open land into built-up land is influenced by lithology.

Factors related to topography and slope gradients have influenced the potential for and development of tourism in the OOW. Tourist attractions have differing morphological conditions and accessibilities. Different landforms have varying hazard potentials due to their distance from the eruption center and slope gradients (Ashari, 2018). The tourism areas in OOW are evenly distributed in the upper, middle and lower parts of the watershed (Figure 9 and Table 10).

It was found that the tourism areas were located at elevations ranging from 2883 meters to 250 meters above sea level, with steep slopes in the upper region, which is part of the volcanic hazard-prone areas of Mount Merapi. The tourism in this region is known as disaster tourism. Based on the disaster prone map (Figure 10), the distribution of tourism areas was concentrated in the volcano hazard-prone areas, mainly clustering in KRB III (Table 11) This area should have been vulnerable for use as a residential and activity zone because it has a high potential for pyroclastic flows, lava flows, volcanic bombs, toxic gases and new collapses. However, due to various beliefs, people have long settled on the slopes of Mount Merapi, and the population in this area continues to grow (Widodo et al., 2024). In the Regent’s Regulation No. 20 of 2011 on the Merapi Volcanic Hazard Zone, KRB III is indeed stated to be unsuitable for residential or tourism development because it could endanger people or visitors during disasters. Unfortunately, some people are still exploiting KRB III as a profitable tourist attraction, while people choose to live in the zone to be close to the source of their livelihoods. Several tourist attractions in KRB III are still facing issues related to permits from the Sleman Regency Government. While tourism near Mount Merapi had boosted the economy, the safety of visitors and the local community needs to be prioritized.

As shown in Figure 10, the change in land-use to built-up areas within the disaster-prone area (KRB) of Mount Merapi in 2034 appears to be very extensive. If this trend continues and encroaches further into the disaster-prone radius, the consequences will be catastrophic, potentially resulting in significant casualties. Therefore, a disaster mitigation-based policy approach is necessary and should be implemented by various relevant stakeholders. One of the key measures is disaster preparedness training and readiness programs for local communities, particularly for tourism operators located in disaster-prone areas. The Regional Disaster Management Agency (BPBD) should intensify such efforts, as they play a critical role in ensuring public safety. Training and evacuation drills for tourism operators who act as first responders during evacuations must be regularly held to ensure that evacuation procedures are conducted in a systematic and effective manner in the event of a disaster.

Additionally, the Regional Development Planning Agency of Sleman Regency should play an active role in disaster mitigation-based spatial planning. This includes restricting development in high-risk zones (red zones), designing emergency evacuation routes, and relocating vital infrastructure to minimize the direct impact of a potential Mount Merapi eruption, particularly in response to the projected increase in built-up areas within tourism zones by 2034. The government could also implement a disaster-mitigation-based tourism concept, whereby visitors not only experience the natural beauty of the area, but also receive education on potential disasters and appropriate response measures. Tourism development should not only prioritize economic aspects, but also ensure safety and environmental sustainability. Such sustainability should be supported by strict regulations regarding land use and tourism facility development to ensure alignment with disaster mitigation strategies. The findings of this study are in line with the research conducted by (Muthohar et al., 2020), who found that some tourist attractions in Mount Merapi volcano KRB III  lack disaster evacuation signs.

However, efforts towards disaster mitigation in terms of accessibility to tourism areas have been initiated, but are hindered by signs and potholes in the roads. Moreover, there is still no well-designed disaster management in the tourism areas. Therefore, there is a need for both structural and non-structural disaster mitigation initiatives in the area, tailored to the morphological characteristics and hazard potentials in order to preserve the terrestrial ecosystem through sustainable land and tourism management, in line with goals 11 and 15 of the Sustainable Development Goals. Consequently, the research findings could serve as a basis for local governments and tourism area managers to formulate land-use policy guidelines integrated with effective regional spatial planning, prioritizing economic growth that synergizes with the considerations of volcanic hazard-prone areas. The study has several limitations that should be considered when interpreting the results. First, it utilizes several Landsat sensors, namely Landsat 7, 8 and 9, each of which has different spectral bandwidths for land use interpretation using natural color composite imagery (red, green and blue).

Therefore, additional processing related to data normalization for each image scene is required when conducting more in-depth analyses involving the spectral values of the imagery, such as the calculation of vegetation density indices (e.g., NDVI). Second, although DEMNAS provides relatively high-resolution national elevation data, micro-topographic variations and rapid post-eruption morphological changes may not be fully captured, particularly in areas affected by newly deposited volcanic materials or intensive erosion processes. Third, the tourism risk analysis in the study focuses primarily on spatial exposure within KRB, but does not quantitatively measure vulnerability components such as the structural resilience of buildings, evacuation preparedness, visitor density, or the adaptive capacity of local communities. In addition, social, cultural and economic factors influencing tourism development decisions are represented only indirectly through spatial variables. Future research could incorporate more dynamic hazard simulations, multi-criteria risk indices, and participatory approaches to enrich sustainable tourism planning analysis in volcanic hazard-prone regions.