Landscape dynamics and its related factors in the Citarum River Basin: a comparison of three algorithms with multivariate analysis

Abstract

Landscape change is intricately linked to natural resource utilization. Landscape dynamics are closely tied to land use and land cover (LULC), serving as a representation of ecosystems and human activities. In the Citarum River Basin, Indonesia, a comprehensive approach is necessary to comprehend landscape dynamics as a manifestation of human interaction with the environment. This research aims to analyze landscape dynamics and its factors that can significantly drive changes. We focused on the Cirasea Watershed, which serves as an upper region of the Citarum River Basin. Data was acquired from Landsat-series imageries from 1993 to 2023, and LULC analyses were conducted using classification and regression trees (CART), random forest (RF), and support vector machine (SVM). We analyzed seven independent variables, including slope (X1), elevation (X2), main river (X3), population (X4), central business district (X5), distance from the past settlements (X6), and accessibility (X7) using multivariate analysis. This research found that RF stands out as the optimal choice for LULC analysis over CART and SVM because it had the highest overall accuracy and overall kappa (0.91–0.92, 0.88–0.89). Notably, there was a substantial 273.43% increase in built-up areas, coupled with significant declines in plantations and horticultures. LULC changes was more pronounced in the lower areas near Bandung City. LR model highlighted X1, X3 and X6 as the significant driving forces for built-up areas expansion (r-square 0.44 with p-value < 0.01 and 95% confidence level). Without effective spatial planning, flat areas near rivers and past settlements have the greatest potential for LULC changes.

RESEARCH HIGHLIGHT

• Landscape dynamics were detected best by random forest (RF), rather than classification and regression trees (CART) and support vector machine (SVM).

• Landscape dynamics were dominated by increasing built-up areas with forests preservation, where plantations and horticultures were significantly reduced.

• Driving factors for landscape dynamics were slope, main river, and distance from the past settlements based on the logistic regression tests.

Keywords:

• Cirasea watershed
• google earth engine
• logistic regression
• LULC
• random forest

Introduction

Landscape dynamics reflected in changes in land use have been proven to have changed environmental quality in many countries (Roy et al. 2022). Landscape dynamics occur massively in many countries pursuing economic growth and development targets (Rodríguez-Rodríguez et al. 2019; Verma and Raghubanshi 2019). Landscape, expressed through land use and land cover (LULC), represents the interaction between environment and human activities. (Kindu et al. 2015). In developing and emerging countries, landscape dynamics typically result in decresing dense vegetation cover, transforming it into human-altered lands, including agricultural, settlements, infrastructure, and industrial areas (Tariq and Mumtaz 2023). At present, landscape management remains highly practical, as regional government-level spatial planning policies are susceptible to modifications driven by economic interests. These changes may promote swift growth in local revenue (Sugito et al. 2023).

Understanding landscape dynamics requires recognizing that changes in LULC are crucial for our future (Dede et al. 2022a). Continuously monitoring the dynamics over an extended period is valuable not just for detecting alterations but also for understanding the spatial and temporal trends of each LULC types. This enables the development of effective measures to prevent any substantial disruption to the environmental quality caused by these changes (Sil et al. 2017; Mello et al. 2020; Daba and You 2022). The rapid development of computing technology has enabled us to make various estimates based on past and current data even today (Widiawaty, Ismail, et al. 2020). In order to achieve sustainable development objectives, a comprehensive understanding of the landscape is essential because various planning regulations related to the development process and environmental management cannot be separated from LULC (Hudalah and Woltjer 2007; Nugroho et al. 2017; Wijaya et al. 2017).

This occurrence underscores the necessity for conducting more comprehensive research on landscape dynamics, with a focus restricted to a watershed or river basin as an ecoregion (Dede, Sunardi, Lam, et al. 2023). Currently, many landscape dynamics studies use formal administrative units, for example, the state down to the village or neighborhood (Widiawaty, Nandi, et al. 2020). This situation can be understood because the development process is divided by regional interests by the local authorities (Rukmana 2015; Setiadi et al. 2023a). Dede et al. (2022b) explored the combination of rural-urban interaction with administrative boundaries in a study focusing on landscape dynamics within ecoregional units. On the other hand, studies on landscape dynamics using watersheds are generally unable to understand the processes and factors that influence changes, as reported by Berget et al. (2021) in Mexico, Gebreslassie (2014) in Ethiopia, Prashanth et al. (2023) in India, and Myers et al. (2023) in the United States. Landscape dynamics, usually characterized by anthropogenic land expansion, is undoubtedly a problem because it able to reduce ecosystem services (Dede, Susiati, Widiawaty, et al. 2023).

Processes in landscape dynamics need to be deeply studied to understand the spatiotemporal factors as drivers, both from biogeophysical and socio-economic aspects (Dede 2020). In Indonesia, several studies focused on watershed landscapes merely discuss their dynamics without delving into the underlying driving forces (Setyorini et al. 2017; Ridwansyah et al. 2020; Djoharam et al. 2022). Many critical watersheds in Indonesia are characterized by floods, landslides, degraded land, and droughts (Susanti et al. 2019). The Citarum River Basin is often in national and international spotlight because this ecoregion is known as the development axis of West Java, Indonesia. Previous studies have focused on landscape changes and have not explored the driving factors in the Citarum River Basin (Sunardi et al. 2022; Ismail et al. 2024). In the other hand, a study from Widiawaty et al. (2019) focused solely on urban sprawl and its typology. Understanding the driving factors that interact with LULC is essential for understanding spatial interactions.

Different from previous studies, our research aims to analyze landscape dynamics and driving factors in the Citarum River Basin. This study will produce a landscape dynamics model as a valuable input to environmental management. A unique approach was complemented by three algorithms: classification and regression trees (CART), random forest (RF), and support vector machine (SVM). These algorithms have different characteristics which require careful consideration of the Citarum River Basin. RF enhances urban vegetation mapping accuracy, surpassing traditional maximum likelihood classifiers, although sensitivity to sampling design exists (Feng et al. 2015; Belgiu and Drăguţ 2016; Georganos et al. 2021). In the other hand, CART classifies and resolves spectral data problems through a tree-structured approach, but bias can arise from data resampling and recursive partitioning (Loh 2008; Ahmadlou et al. 2022). Furthermore, SVM addresses image classification fragmentation and low accuracy issues (Jia and Wang 2011; Nandam and Patel 2022), yet an increase in features, especially with a small training sample, may significantly reduce classification accuracy (Pal and Foody 2010). The chosen algorithm, coupled with multivariate investigations, was employed to reveal driving factors using logistic regression (LR) and geographically weighted logistic regression (GWLR). This method not only enhances creative contributions to geospatial research but also holds specific significance in the Citarum River Basin management.

Methodology

This study was conducted in a sub-region of the Citarum River Basin, specifically in the Upper Citarum section, West Java Province. Upper Citarum serves as a crucial water catchment area with a strategic role in the development of West Java, particularly in Greater Bandung (Khairunnisa et al. 2020; Fitriana et al. 2023). This region is divided into smaller ecoregions, such as Cikapundung, Cirasea (Cisangkuang), Cihaur, Ciwidey, Cisangkuy, Ciminyak, and Citarik segments, with their respective watersheds (Marko et al. 2021). Our research focuses on the Cirasea or Cisangkuang Watershed, as shown in Figure 1. Cirasea Watershed covers an area of 366.67 km², and the main river flows from the mountainous area in the south. This area was once the focus of the Citarum Harum program (Handayani 2019) and is known as the textile industries center in Greater Bandung. In addition, Cirasea Watershed has developed rapidly as a development pole. Therefore, inevitable landscape changes include the expansion of residential areas and infrastructure supporting Majalaya’s industrial activities (Fathunnisa and Salami 2015; Sunardi et al. 2021; Burnama et al. 2022).

Figure 1. Cirasea Watershed as the research focus with field survey locations.

Data acquisition

Our study primarily examined landscape dynamics, with a focus on the classification of LULC into seven distinct types: water bodies, empty land (non-vegetated), agriculture, built-up area, forests, shrubs, grassland, plantations and horticultures (Widiawaty, Ismail, et al. 2020). Table 1 presents more details regarding the research variables and indicators. LULC data refers to the interpretation of Landsat-series satellite imageries originating from the United States Geological Survey (USGS) for a 10-year period from 1993 to 2023. This interval was selected following the period for clear visible changes (Dede, Susiati, Widiawaty, et al. 2023). This research contains data originating from the dry season, which is free from cloud cover, precisely from April 1st to October 30th each year (Susiati, Dede, et al. 2022; Oktavia et al. 2023; Rachmadita et al. 2023;). The information underwent remote sensing processing to generate eight LULC types. Satellite imagery was interpreted using a supervised classification method (Widiawaty 2019). CART, RF, and SVM were selected as the algorithms for classifying remote sensing images (Lam and Noor 2018; Fikri et al. 2021; Basheer et al. 2022). In this context, CART pertained to 10 lead nodes, RF involved 50 trees, and SVM utilized linear kernel that considers regional characteristics.

Table 1. Variables to dataset information.

Variable	Parameter	Information
Landscape dynamics	LULC types	Water bodies, non-vegetated, agriculture, built-up area, forests, shrubs, grasslands, as well as plantations and horticultures
	Changes for each LULC	Changes in area and percentage (%) every 10 years from 1993 to 2023
Driving factors for landscape changes	Biogeophysical factors	Slope (X1), elevation (X2) and main river (X3)
	Socio-economic factors	Population (X4), central business district (X5), distance from the past settlements (X6), and accessibility (X7)

These classifications use 80 regions of interest (ROI) for each LULC type with 80% training and 20% testing. The selection of training data is based on high-resolution satellite imagery from PlanetScope and the researcher’s field experience from the previous year, ensuring purposive selection of homogeneous regions (pixels) (Patnala et al. 2023). In general, the training is larger than the testing to train algorithms to accurately detect LULC (Feizizadeh et al. 2023). All bands from three satellite images from Landsat-5 TM (1990), Landsat-7 ETM+ (2003 and 2013), and Landsat-8 OLI (2023) are used. LULC data was processed using Google Earth Engine as a cloud-based platform for processing satellite imagery and QGIS version 3.28 (Firenze) after the image selection stages (cloud cover less than 10%), cloud masking, and preprocessing to display earth surface phenomena with minimal errors (Sano et al. 2007; Zhu and Woodcock 2014; Madzen and Lam 2017). In addition, independent variable data, which act as driving forces, comes from various credible sources. To facilitate further analysis, geospatial data must be used (Pebesma and Bivand 2023). The slope comes from height difference information, which is analyzed using the slope analysis method (Böhme et al. 2011), whereas other data must be analyzed using buffering and density analysis (Susiati, Widiawaty, et al. 2022). Table 2 presents details of data processing on biogeophysical and socioeconomic factors in order to provide information on driving factors. The independent variables were validated in prior landscape dynamics studies in Indonesia and served as significant geospatial driving factors (Sodikin et al. 2017; Widiawaty, Ismail, et al. 2020; Dede et al. 2022b).

Table 2. Data acquisition for biogeophysical and socio-economic factors.

Factor	Method	Information	Data source
X1	Slope analysis	Comparison between flat distance and height difference at two different coordinate points (%)	Geospatial Information Agency (BIG)
X2	Physiographic	The height above average sea level – meters above sea level (masl)	Geospatial Information Agency (BIG)
X3	Buffering	Distance data of a place from the main river, in meters (m)	Indonesian Ministry of Public Works and Housing of Indonesia
X4	Density analysis	Population density in a particular area refers to population per square kilometers (km^2)	Indonesian Statistical Agency (BPS); Indonesian Ministry of Home Affairs; WorldPop-University of Southampton
X5	Buffering	Distance data from the central business district that the government has determined, represented in meters (m)	Regional Development Planning Agency of Bandung and Garut
X6	Buffering	Flat distance from the existing settlement matrix and patch in the previous period, represented in meters (m)	LULC Maps from satellite imageries
X7	Road density analysis	The comparison between the length of roads (arterial, collector and local) in a particular area, represented in units of km per km^2	Geospatial Information Agency (BIG); HOT Open Street Maps; Ministry of Public Works and Housing of Indonesia

Data validation

To ensure the quality of landscape information, validation of the interpretation results involves 86 ground-check points (field survey) adjusted to their allocation value in the Cirasea Watershed in 2023 for each LULC type (Figure 2). The ground-check points were selected through a purposive sampling method, considering the distribution and accessibility of locations. Essential information includes geographic coordinates, photographs (sketches), and timestamps for data collection (Murayama 2012; Chaminé et al. 2021). In the validation phase, both overall accuracy and overall kappa were employed. Overall kappa demonstrated superior discernment as a measure for accurate classification in remote sensing data, exhibiting enhanced discrimination between classes compared to overall accuracy. (Li et al. 2019). Landscape information was declared valid after testing the overall accuracy and overall kappa with thresholds between 0.70 and 0.96 with a confidence level of 95% (Prasad et al. 2022). Determining the threshold value maintains ideal landscape data and avoids underfitting or overfitting (Montesinos-López et al. 2022; Widiawaty et al. 2022). The series of landscape data acquisition and validation is presented in Figure 2.

Figure 2. Framework to obtain LULC data as a landscape representation.

Data analysis

Landscape dynamics were analyzed by comparing LULC information over 10 years following the gain and loss method (Widiawaty 2019). These landscape changes are known from the changes in each LULC area. An increase trend signifies a gain, whereas a decrease indicates a loss. Analyzing these gains and losses provides change values in area units (km²) or percentages (%) (Ismail et al. 2020). This descriptive landscape change pattern can be clearly depicted using a trend line. Changes in each LULC were represented in binary form, categorizing them as either changed or unchanged. These binary representations served as the dependent variable (Toma et al. 2023). The effect was assessed based on the coefficient of determination and the p-value at a 95% confidence level (Danapriatna et al. 2023). To prevent symptoms of multicollinearity, typically indicated by an r-value exceeding 0.90 with statistical significance, it is crucial for the various factors serving as independent variables to establish their relationships with one another (Widiawaty et al. 2022).

The influence of these factors on landscape changes in the Cirasea Watershed was determined by LR and GWLR analysis using PSPP version 1.6.2. LR models the relationship between independent and dependent variables using categorical or binary data (Tsangaratos and Ilia 2016). GWLR involves geographic factors that combine non-stationary and categorical parameters (Caraka and Yasin 2017; Nkeki and Asikhia 2019). In this investigation, the utilization of binary numbers derived from LULC change data is well-suited for analysis employing the regression method outlined in Equation 1(1) $$\mathrm{\pi }\mathrm{ }\left({\mathrm{X}}_{\mathrm{i}}\right)=\frac{\mathrm{\text{exp}}\hspace{0.17em}\left({\sum }_{\mathrm{j}=0}^{\mathrm{p}}{\mathrm{\beta }}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}})\right)}{1+\mathrm{\text{exp}}\hspace{0.17em}\left({\sum }_{\mathrm{j}=0}^{\mathrm{p}}{\mathrm{\beta }}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}){\mathrm{x}}_{\mathrm{\text{ji}}}\right)}$$(1) . The GWLR model employed in this study incorporates both adaptive-Gaussian and fixed-Gaussian kernels. The model is deemed superior if the values for deviance (Dev), degrees of freedom (DoF), Dev per DoF, and Akaike information criterion (AIC) are lower compared to the LR model without geographic weighting as per Equation 2(2) $$\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{\text{1‐}{\mathrm{P}}_{\mathrm{i}}}\right)\text{=}\mathrm{\alpha }+{\mathrm{\beta }}_1\mathrm{ }{\mathrm{X}}_1+{\mathrm{\beta }}_2\mathrm{ }{\mathrm{X}}_2+\dots .{\mathrm{\beta }}_{\mathrm{n}}\mathrm{ }{\mathrm{X}}_{\mathrm{n}}\mathrm{ }$$(2) (De and Acquah 2010; Dede et al. 2022a). (1) $$\mathrm{\pi }\mathrm{ }\left({\mathrm{X}}_{\mathrm{i}}\right)=\frac{\mathrm{\text{exp}}\hspace{0.17em}\left({\sum }_{\mathrm{j}=0}^{\mathrm{p}}{\mathrm{\beta }}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}})\right)}{1+\mathrm{\text{exp}}\hspace{0.17em}\left({\sum }_{\mathrm{j}=0}^{\mathrm{p}}{\mathrm{\beta }}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}){\mathrm{x}}_{\mathrm{\text{ji}}}\right)}$$(1) where X_ji is the observed value of the predictor variable at location (ui, vi), βj (ui, vi) is the regression coefficient for each location (ui, vi), and p is the number of predictor variable parameters. (2) $$\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{\text{1‐}{\mathrm{P}}_{\mathrm{i}}}\right)\text{=}\mathrm{\alpha }+{\mathrm{\beta }}_1\mathrm{ }{\mathrm{X}}_1+{\mathrm{\beta }}_2\mathrm{ }{\mathrm{X}}_2+\dots .{\mathrm{\beta }}_{\mathrm{n}}\mathrm{ }{\mathrm{X}}_{\mathrm{n}}\mathrm{ }$$(2) where ln is the natural logarithm, p the probability value, α the regression constant, and β the predictor coefficient of the independent variable.

The GWLR and LR models were then tested for feasibility using the −2 log-likelihood value, Nagelkerke test, and Hosmer-Lemeshow test provided by statistical computing software (Park 2013; Smith and McKenna 2013). The 2 log-likelihood value and the Nagelkerke test (r-square) are declared significant if the p-value 0.05, whereas the Hosmer-Lemeshow test states otherwise (Widarjono 2015; Barbić and Palić 2016; Ahmed 2017). These regression models must meet the heterogeneity aspect by referring to the Breusch-Pagan test and avoid multicollinearity using the variance inflation factor (VIF) method – it must be no more than 10.00 (Vu et al. 2020; Solekha and Qudratullah 2022). The research variables showed heterogeneity because the Breusch-Pagan value reached 1257.71 (p-value 0.00 < 0.05) and was free from multicollinearity. The VIF reached 4.30, making it suitable for further analysis.

Results

Three supervised classification algorithms yielded varying results regarding LULC in the Cirasea Watershed. RF exhibited the highest level of accuracy in detecting the landscape, as delineated in Table 3. RF demonstrated superiority over both CART and SVM, as evidenced by the utilization of ‘testing data’ in Google Earth Engine and ‘field survey data’ conducted in 2023, achieving an overall accuracy of 0.91–0.92 and an overall kappa of 0.88–0.89, respectively. SVM exhibited a lower level of accuracy in LULC detection compared to alternative methods due to its tendency to over-detect forests and misclassify rice fields as plantations and horticultures (Figure 3). Conversely, SVM also demonstrated reduced precision when identifying grassland, marsh, and non-vegetated areas. Meanwhile, CART exhibited a spatial detection pattern for LULC with relatively favorable outcomes, achieving accuracy levels between those of RF and SVM. However, CART displayed a narrower tendency towards over-detection in grassland, marsh, and non-vegetated areas. CART was not the least effective in detecting LULC associated with human cultivation, particularly on agricultural lands, built-up areas, plantations, and horticulture. However, initial observations suggest that CART yields results akin to those of RF. Discrepancies in detection among grassland, marsh, and non-vegetated areas stem from the closely clustered reflectance values of these entities, rendering them difficult to distinguish without field verification (Tong et al. 2018; Gao et al. 2023).

Figure 3. Comparison of LULC results of three classification algorithms.

Table 3. Accuracy results between CART, RF and SVM.

Classification method	Google Earth Engine		Field survey
	Overall accuracy	Overall Kappa	Overall accuracy	Overall Kappa
CART	0.82	0.77	0.76	0.68
RF	0.92	0.89	0.91	0.88
SVM	0.82	0.76	0.65	0.55

Consistency was observed in the results concerning the identification of developed regions within the Cirasea Watershed. Patterns, distributions, and extents exhibited no significant disparities, suggesting that various supervised classification algorithms adeptly detect landscapes heavily influenced by human activity. RF emerged as the most suitable algorithm for detecting landscape dynamics in this watershed, considering the minimal thresholds for overall accuracy and kappa accuracy. RF effectively discerns LULC features such as agricultural land, built-up areas, forests, as well as plantations and horticultures, which represent the predominant landscapes therein. RF demonstrated adaptability in detecting landscapes ranging from highly natural to those dominated by human activity. Discrepancies in satellite image values are effectively captured by RF, facilitating the production of precise LULC information (Belgiu and Drăguţ 2016; Santos et al. 2019).

LULC changes

The Cirasea watershed is situated in the southern part of Citarum River Basin. Three decades ago, the landscape primarily comprised agricultural land and forests associated with plantations and holticultures. This landscape configuration has persisted for the past decade (Figure 3); however, by 2023, notable changes have ensued due to the continued expansion of built-up areas, resulting in a drastic alteration exceeding 13% (Table 4). The proliferation of built-up areas has led to the establishment of new villages and housing complexes. Furthermore, it is pertinent to highlight the influence of the textile industry on shaping landscape transformations in the Cirasea Watershed. Built-up areas for residential, amenity, and economic purposes exhibit spatial interconnectedness, with their locations relatively proximate to each other (Widiawaty, Ismail, et al. 2020). The downstream region of the Cirasea Watershed serves as the epicenter for the textile industry’s development in Greater Bandung, evidenced by the presence of multinational manufacturers situated in close proximity to highways and settlements spanning from the eastern (Majalaya) to western (Baleendah) sides (Handayani 2019; Ripaldi 2023). This phenomenon has dynamically altered the proportions of the four LULC types during the period spanning 1993 to 2023 (Figure 4).

Figure 4. Distribution of LULC in the Cirasea Watershed from 1993 to 2023.

Table 4. LULC dynamics in the Cirasea Watershed.

LULC	1993		2003		2013		2023
	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%
Agricultural land	77.56	21.15	62.71	17.10	62.65	17.09	70.15	19.13
Built-up areas	12.80	3.49	18.58	5.07	23.00	6.27	47.80	13.03
Forest	71.86	19.60	67.20	18.33	76.31	20.81	88.72	24.19
Grassland	16.64	4.54	4.36	1.19	3.05	0.83	14.15	3.86
Marsh	18.23	4.97	26.99	7.36	25.28	6.90	12.31	3.36
Non-vegetated	17.60	4.80	29.85	8.14	18.55	5.06	7.87	2.15
Plantations and horticultures	151.13	41.21	156.66	42.72	157.45	42.94	124.57	33.97
Water bodies	0.88	0.24	0.35	0.09	0.40	0.11	1.14	0.31
Total	366.69	100.00	366.69	100.00	366.69	100.00	366.69	100.00

Gain and loss analysis reveals a more than 270% increase in built-up area over the past three decades, accompanied by reductions in agricultural land, grassland, marsh, non-vegetated land, plantations, and horticulture (Table 5). Notably, forests have witnessed an increase in both area and proportion within the Citarum River Basin since 2003, attributable to landscape restoration initiatives and structured governmental efforts to mitigate critical land degradation (Pitaloka et al. 2020; Ayyasy et al. 2021). This trend contrasts with agricultural land, which experienced significant changes during 1993–2003, primarily characterized by reductions to accommodate various built-up areas. Over the past decade, the built-up area has expanded by over 100%, underscoring the heightened demand for land for residential and industrial purposes. Additionally, water bodies have expanded, facilitated by the construction of irrigation canals and small dams to cater to the requirements of industrial and agricultural activities for raw water.

Table 5. Gain and loss of LULC in 1993–2023.

LULC	1993–2003		2003–2013		2013–2023		1993–2023
	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%
Agricultural land	−14.85	−19.15	−0.06	−0.09	7.50	11.96	−7.41	−9.56
Built-up areas	5.78	45.14	4.42	23.82	24.79	107.80	35.00	273.43
Forest	−4.66	−6.48	9.10	13.55	12.41	16.26	16.86	23.46
Grassland	−12.28	−73.79	−1.31	−29.99	11.10	363.54	−2.48	−14.93
Marsh	8.76	48.07	−1.71	−6.33	−12.98	−51.33	−5.92	−32.49
Non-vegetated	12.25	69.58	−11.30	−37.87	−10.68	−57.58	−9.73	−55.31
Plantations and horticultures	5.53	3.66	0.79	0.50	−32.88	−20.88	−26.56	−17.57
Water bodies	−0.54	−60.79	0.05	15.81	0.74	185.06	0.26	29.43

Driving factors in landscape dynamics

The built-up area in the Cirasea Watershed continues to expand, nearly tripling between 1993 and 2023. Investigating this phenomenon using the LR and GWLR models to discern the spatial-temporal factors influencing its dynamics is of considerable interest. The driving forces are categorized into three biogeophysical and four socio-economic factors with distinct spatial distribution characteristics, as depicted in Figure 5. As independent variables, these driving factors can construct a model with an accuracy level of up to 94.70%. This figure is significant as it yields −2 log-likelihood value lower than that of the model without log likelihood, with a difference of 222.70 (Table 6). Among the three logistic models examined, all produced significant r-square values ranging from 0.44 to 0.61 (p-value 0.01, 95%), indicating that the independent variables can elucidate landscape dynamics, particularly regarding open-land expansion in the Cirasea Watershed (Table 7). The Hosmer-Lemeshow test yielded results of 4.94 (p-value > 0.76, 95%) for both GWLR models in explicating the driving forces behind LULC dynamics (Table 8). However, these two GWLR models cannot be employed for LULC analysis due to the absence of spatial variability as a local geographical variation (Mayfield et al. 2018).

Figure 5. Biogeophysics and socio-economics drivers of LULC changes.

Table 6. Validation of the LR model.

Model stage	Overall accuracy (%)	−2 log-likelihood
Without independent variables	94.90	567.88
Using independent variables	94.70	345.18

Table 7. Hosmer-Lemeshow Test and simultaneous determination of driving forces.

Model	Chi-square	dF	Nagelkerke r-square	p-Value	Epsilon
LR	4.94 (0.76)	8.00	0.44	<0.00	0.56
GWLR (Adaptive Gaussian)			0.61	<0.00	0.39
GWLR (Fixed Gaussian)			0.52	<0.00	0.48

Table 8. Comparison of various logistic regression models and their variances.

Parameter	GWLR		LR
	Adaptive Gaussian	Fixed Gaussian
Bandwidth size	6.86	133.00	0.00
Dev	226.25	243.38	265.02
DoF	1,368.11	1,363.22	1,381.00
Dev per DoF	0.16	0.18	0.19
AIC	268.69	269.42	281.12

Note: positive value of AIC suggests no spatial variability.

The positive AIC values, with magnitudes almost identical between spatial and non-spatial regression models, imply that classic LR is better suited to elucidate this phenomenon (Nakaya 2012; Mondal et al. 2015). The high r-square values (0.52 and 0.61) in the logistic regression models and the two GWLRs raise concerns regarding overfitting. In logistic regression models, low r-square values are common, typically around 0.40 if the model meets goodness of fit (Hosmer et al. 2013). The lack of spatial variability stems from the extensive and intensive growth of built land in the downstream part of the Cirasea Watershed, driven by numerous forces. Consequently, the LR model emerges as the most suitable model, with a driving force determination of 44% and omnibus test results of 222,697 (p-value 0.00 and dF 8.00, 95%). According to the LR model, slope (X1), the main river (X3), and distance from past settlements (X6) have emerged as significant driving forces for built-up area expansion in the Cirasea Watershed. New built-up areas tend to appear in areas with gentle slopes near surface water sources (rivers), characterized by a negative coefficient and close proximity to settlement patches and matrices from the previous period. The LR model for landscape dynamics is expressed by Equation 3(3) $$\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{\text{1‐}{\mathrm{P}}_{\mathrm{i}}}\right)\text{= 0.5151}-0.0615\mathrm{ }{\mathrm{X}}_1+{0.0005\mathrm{ }\mathrm{X}}_2-0.0009\mathrm{ }{\mathrm{X}}_3-0.0001\mathrm{ }{\mathrm{X}}_4-0.0001\mathrm{ }{\mathrm{X}}_5-0.0052\mathrm{ }{\mathrm{X}}_6+{0.1472\mathrm{ }\mathrm{X}}_7$$(3) , as derived from Table 9.

Table 9. Parameters of the LR model.

Variable	B	SE	Wald	p-Value	Exp(β)
X1	−0.0615	0.02	6.88	0.01	0.94
X2	0.0005	0.00	0.05	0.82	1.00
X3	−0.0009	0.00	8.19	0.00	1.00
X4	−0.0001	0.00	1.56	0.21	1.00
X5	−0.0001	0.00	0.66	0.42	1.00
X6	−0.0053	0.00	26.89	0.00	0.99
X7	0.1472	0.19	0.58	0.45	1.16
Constant	0.5151	1.24	0.17	0.68	1.67

Note: β = independent variable predictor coefficient, SE = standard errors of the regression, and Exp = exponent value of β or odds ratio.

Within the LR model, only three driving forces exhibit significance, with no significant correlation observed with an r-value less than 0.10 (p-value > 0.05, 95%), as illustrated in Table 10. The correlation between slope (X1) and the main river (X3) is very weak (0.08), as is the correlation between slope (X1) and the distance from past settlements (X6) (0.03), and between the main river (X3) and the distance from past settlements (X6) (0.02). Overall, the correlation pattern between independent variables is notably low, with most values falling below 0.40, further indicating the absence of multicollinearity (Das 2019; Frost 2020; Patil and Franken 2021). Correlations ranging from weak to high among independent variables that were not significant for the LR model are observed. The highest correlation is between the distance from the central business district (CBD) (X5) and elevation (X2), and population (X4), reaching 0.83 and 0.58, respectively.

Table 10. Correlation between driving forces as independent variables.

Variable	X1	X2	X3	X4	X5	X6	X7
X1	1.00	–	–	–	–	–	–
X2	−0.13	1.00	–	–	–	–	–
X3	0.08	−0.17	1.00	–	–	–	–
X4	0.04	−0.37	0.11	1.00	–	–	–
X5	0.01	−0.83	0.04	0.58	1.00	–	–
X6	0.02	−0.01	0.02	0.27	−0.01	1.00	–
X7	0.11	0.38	0.07	−0.28	−0.42	0.06	1.00

(3) $$\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{\text{1‐}{\mathrm{P}}_{\mathrm{i}}}\right)\text{= 0.5151}-0.0615\mathrm{ }{\mathrm{X}}_1+{0.0005\mathrm{ }\mathrm{X}}_2-0.0009\mathrm{ }{\mathrm{X}}_3-0.0001\mathrm{ }{\mathrm{X}}_4-0.0001\mathrm{ }{\mathrm{X}}_5-0.0052\mathrm{ }{\mathrm{X}}_6+{0.1472\mathrm{ }\mathrm{X}}_7$$(3)

Discussion

The dynamics of Land Use and Land Cover (LULC) exhibit three primary patterns: increasing, decreasing, or stagnant (Figure 6), contingent upon changes in area over the past 30 years. Agricultural land has dwindled by several dozen km², yet in the last decade, a reversal has occurred, resulting in a slight increase, though not as significant as in 1993. This trend is also noticeable in grasslands, albeit with a narrower area. Marsh and non-vegetated areas underwent a pattern characterized by an initial upward trajectory followed by a sharp decline. This shift stemmed from the repurposing of these lands for other uses, previously considered abandoned (Jalayer et al. 2022). Over the initial two decades, plantations and horticultures experienced growth, but in the last decade, there has been a notable loss of land, exceeding 20 km². The dynamics of LULC in the Cirasea Watershed indicate an increasing trend in built-up areas. This escalation is commonplace, reflecting similar phenomena observed elsewhere, particularly in peri-urban regions (Dede et al. 2019; Admasu et al. 2023a, 2023b).

Figure 6. The trend of LULC changes in Cirasea Watershed.

Over the past decade, built-up areas have expanded by over 20 km², with exponential growth patterns attributed to industrialization and real estate development, particularly to accommodate workers originating from outside the region (Parveen et al. 2022). Ultimately, these workers chose to settle near industrial sites to mitigate commute times, as the main road from Baleendah to Majalaya often experienced traffic congestion during peak commuting hours, notwithstanding occasional flooding (Junnaedhi et al. 2017; Setiadi et al. 2023b). Conversely, the expansion of built-up areas follows a pattern not solely reliant on forests or agricultural land but can also arise from plantations and horticultures, and may even encroach upon previously marshy or non-vegetated areas (Figure 7). The sustainability of these two land types is relatively assured owing to forest conservation efforts and enhanced irrigation practices. Forest protection is rigorously enforced by both governmental and community-led initiatives through community-based management (Idris et al. 2019). Since 2003, forests have continued to expand, surpassing the area recorded in 1993, a period marked by extensive deforestation in Indonesia and other developing countries due to the extensive utilization of tropical forest resources (Liping et al. 2018)

Figure 7. Landscape dynamics model in the Cirasea Watershed with its significant driving forces and potential impact.

LULC in the Cirasea Watershed will continue to be dynamic as because the population of working-age is still to grow. There will be a definite need for built-up area to accommodate economic activities and housing (Lam et al. 2020). The peri-urban landscape pattern encourages the emergence of built-up area with a sprawl pattern (Widiawaty et al. 2019), especially the Cirasea Watershed right to the south of Bandung City, which is connected via arterial and collector roads. Cirasea’s position as part of Megapolitan Bandung Raya while maintaining forests cannot be separated from the basin-like physiographic ‘Cekungan Bandung’ (Mulyadi et al. 2020). The expansion of built-up area in Cirasea is predominantly located in flat to undulating areas, the hilly and mountainous areas at the top of the Bandung Basin remain forested due to spatial zoning (Truax et al. 2015; Aisharya et al. 2022). Expansion of built-up areas tends to be close to the main rivers, which raises concerns about this. The Cirasea River has the potential to become a dumping ground for household and industrial waste, ultimately disrupting the aquatic ecosystem of the Cirasea Watershed (Alamanda et al. 2016).

In the forthcoming years or decades, future research is imperative to examine the potential landscape dynamics of the Cirasea Watershed and the Citarum River Basin. A profound understanding of various driving factors forms the essential foundation for designing a dynamic model capable of accommodating scenarios tailored to current conditions, such as socio-economic projections and government spatial planning directions (He et al. 2008). Furthermore, the scope needs expansion by incorporating new independent variables such as land prices, vulnerability to various disasters, and public perceptions of space (Widiawaty and Dede 2018). To enhance the richness of LULC types, landscape dynamics analysis requires more detailed data sources, including the utilization of unmanned aerial vehicles (UAVs) (Gomez et al. 2022). Exploring optimized and integrated LULC algorithms can be a crucial step in refining landscape classification methodologies (Escobar-López et al. 2024). Overall, this research marks the beginning of further development, and the current limitations, which have not reached a very detailed scale, present opportunities for further studies. Additionally, future research should explore the implications of landscape dynamics on environmental parameters, encompassing air, water, soil, and other species, transcending their impact solely on humans (Werner and McNamara 2007; Grimaldi et al. 2014; Melesse and Abtew 2016; Olson et al. 2017).

Conclusion

The RF algorithm emerges as the most suitable for analyzing landscape dynamics in the Cirasea watershed. Overall accuracy and kappa values are higher for RF than for CART and SVM. Over the 30-year period (1993–2023), LULC changes were predominantly characterized by the expansion of built-up areas and a reduction in plantations and horticultures. Simultaneously, the LR model identified all three independent variables as driving forces: slope (X1), main river (X3), and distance from past settlements (X6), influencing 44% of the LULC changes. These variables exhibit negative coefficient notation, signifying that built-up areas tend to manifest in locations that are flat and close to water sources. The increasing forest area in the Cirasea Watershed is encouraging despite pollution resulting from industrial activities and critical lands. However, this study is limited as it examines landscape change phenomena without developing a spatial-temporal prediction model. Efforts to control the expansion of built-up areas must secure structural and non-structural support involving active participation from all stakeholders, including government, society, academics, researchers, industry, and non-governmental organizations. The Indonesian government should broaden the Citarum Harum program to encompass landscaping beyond its previous exclusive focus on the restoration of the river water ecosystem. Local authorities must enforce stringent legal measures for housing and industrial developers, including authorization revocation, imposing high sanctions, and promoting willingness to pay for ecosystem services.

Disclosure statement

No potential conflict of interest was reported by the author(s).