The study encompassed the whole island of Java which is located on the southern part of the Indonesian archipelago (Figure 1). The island spans 132,792 km², extending from 5°52′S, 105°04′E in the northwest to 8°47′S, 114°36′E in the southeast. Administratively, it is divided into six provinces : Banten, DKI Jakarta, West Java, Central Java, Yogyakarta, and East Java. Java has a tropical climate characterised by consistently warm temperatures and high humidity throughout the year, with seasonal variations in rainfall, with most regions receivinge between 2,000 and 4,000 mm of annual precipitation. The wet season typically occurs from November to March, bringing substantial rainfall to much of the island, while the eastern parts of Java tend to remain relatively dry from May to September (Sosilawati et al., 2017). Almost all lowlands are used for agricultural and settlements, while the highlands are dominated by an increasing number of plantation crops (Setiawan & Kunihiko, 2020). The human population continues to grow, along with the expansion of residential areas into non-urban areas (Pravitasari et al., 2024). However, Java still aprovides suitable habitat for six threatened species based on data from the Indonesia citizen science project (Squires et al., 2021).

Data were collected for 13 kingfisher species (the blue-eared kingfisher (Alcedo meninting), the Javan blue-banded kingfisher (Alcedo euryzona), the cerulean kingfisher (Alcedo coerulescens), the rufous-backed kingfisher  (Ceyx rufidorsa), the banded kingfisher (Lacedo pulchela), the Javan kingfisher (Halcyon cyanoventris), the collared kingfisher (Todirhampus chloris), the white-throated kingfisher (Halcyon smyrnensis), the sacred kingfisher (Todiramphus sanctus), the ruddy kingfisher (Halcyon coromanda), the common kingfisher (Alcedo atthis), the stork-billed kingfisher (Pelargopsis capensis) and the black-capped kingfisher (Halcyon pileate)) from the period 2016 to 2023. We combined data from multiple citizen science (CS) platforms, thus obtaining a more balanced and comprehensive perspective (Silva et al., 2023). The CS data expanded survey coverage into remote regions with limited road access, thereby reducing spatial bias (Squires et al., 2021). Data were obtained from the Burungnesia CS project (https://birdpacker.org/burungnesia-kupunesia/); eBird (https://ebird.org/, accessed on 13 March 2024); and the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/, accessed on 04 October 2023). These CS platforms follow standardised protocols for data recording, including the use of complete checklists (Callaghan et al., 2021; Kelling et al., 2019) with data subsequently verified by a team of local experts (Sullivan et al., 2009).

Data filtering began by downloading occurrence data based on human observation or protocol type (traveling, stationary and incidental) (Johnston et al., 2021; Kelling et al., 2019). Occurrence records for kingfisher species were checked for duplicate entries, which were subsequently removed (Chowdhury et al., 2023; Galván et al., 2022) by using the “Delete Duplicate” feature in ArcGIS Pro 3.2. Spatial filtering was then applied in order to retain only one presence point per species within its estimated home range (Lajeunesse & Fourcade, 2023). To reflect biological realism, we considered known variations in home range size among Alcedinidae subfamilies. For example, species in the Halcyoninae subfamily are generally associated with larger home ranges than those in the Alcedininae family, largely due to differences in body size and behavior (Kesler & Haig, 2007a; Musseau et al., 2021). Accordingly, we used one occurrence per grid cell of 7.29 km² for Halcyoninae and 5.74 km² for Alcedinidae (Musseau et al., 2021). We applied a spatial rarefaction filter using home range estimates for Alcedininae as a baseline (Lajeunesse & Fourcade, 2023).

Thirteen kingfisher species have been recorded on Java, of which nine were retained for modeling after filtering: the blue-eared kingfisher (N = 715), the Javan blue-banded kingfisher (N = 38), the cerulean kingfisher (N = 1103), the rufous-backed kingfisher  (N = 144), the banded kingfisher (N = 116), the Javan kingfisher (N = 2946), the collared kingfisher (N = 3671), the white-throated kingfisher (N = 33), and the sacred kingfisher (N = 197). The black-capped kingfisher was excluded due to the availability of only one presence point following filtering, which was insufficient for model training and validation. The ruddy kingfisher, the common kingfisher, and the stork-billed kingfisher were also excluded due to the caution needed when interpreting results from small samples (N < 30) (Breiner et al., 2015; Cao et al., 2024).

To understand the distribution patterns of various kingfisher species in relation to geographic and environmental factors, we selected key environmental variables based on their ecological relevance and common use in habitat assessments (Fang et al., 2025). Those considered for the study included vegetation density, temperature, elevation, vegetation strata, distance to water bodies, proximity to built-up areas, and annual rainfall (Table 1). These data were reprojected and all environment variables were standardised to the same coordinate system (WGS 1984 UTM Zone 49S) with a spatial resolution of 100 x 100 m, and saved in GeoTIFF file format, as required by the random forest model, in order to facilitate interpretation of the model outputs. Prior to model development, multicollinearity among the selected variables was analysed using Pearson’s correlation analysis. A dendrogram derived from hierarchical clustering of correlation distances was used to visualize the structure of the relationship between predictors and to identify clusters of highly-correlated variables. In cases where two or more variables were highly correlated cluster (|r| > 0.70), we retained the variable that was more ecologically meaningful (Achour & Kalboussi, 2020) in relation to the known habitat preferences and requirements of kingfisher species.

Webuilt habitat suitability models using a random forest (RF) algorithm implemented in the randomForest package in R version 4.1.1 (2021-08-10). RF is relatively robust to multicollinearity and overfitting due to its ensemble structure particularly in relation to spatial citizen science data (Imron et al., 2025). The models were built with a collection of decision trees to obtain accurate and stable predictive information . Model robustness was further assessed using repeated cross-validation (Ramampiandra et al., 2023) and multiple evaluation metrics (Figure 2).

Absence data were approximated through pseudo-absence generation, simulating non-detection areas within the environmental space of the study region (Broussin et al., 2024). The RF model, applied to both presence and pseudo-absence samples, produced reliable predictions, thereby maximising the value of the data collected through citizen science initiatives (Bird et al., 2014; O’Neill et al., 2023). The number of pseudo-absence sample data for each kingfisher species was the same as the presence sample points . The approach taken restricts the estimated environmental niche of a species because the model becomes dominated by absence data (Broussin et al., 2024). Pseudo-absence points were randomly generated outside the buffer area using the 'Generate Random Points' tool in ArcGIS Pro. The buffer area was an average dispersal distance of 849 meters (as the radius of the circle), reflecting the typical movement range of kingfisher, which aggressively defend their territories during the breeding season . All background locations were considered to represent the possibility of not finding the species at a particular location and were treated as pseudo-absence data (Sillero et al., 2021; Zhang et al., 2019).

Model performance was evaluated using repeated 5-fold cross-validation implemented through the caret R-package. Such cross-validation is essential for obtaining an estimate of predictive performance and to detect potential overfitting, ensuring that the RF model captures ecologically meaningful relationships rather than noise in the data (Ramampiandra et al., 2023). The discrimination ability of the model was quantified using the area under the curve (AUC). Cross-validation AUC values were reported as mean standard deviation (SD) to reflect model stability across folds. The RF model (ntree = 10,000) was trained on the full dataset to produce continuous habitat suitability prediction. A higher number of trees is beneficial in reducing model variance and ensuring reliable estimation of variable importance . RF performance has been shown to be relatively robust to parameter settings (Sillero et al., 2021). The number of predictor variables considered at each node split (mtry) followed the default RF classification setting, corresponding to the square root of the total number of predictor variables (Albertini et al., 2024) in the randomForest package. Binary presence-absence maps were derived using a data-driven threshold based on receiver operating characteristic (ROC) (Brownscombe et al., 2021). Sensitivity and specificity (TPR + TNR) were calculated for possible thresholds derived from RF predictions. Data-based thresholds, determined by maximising the sum of TPR and TNR, provide objective criteria for converting continuous habitat suitability predictions into biner presence-absence maps (Sahragard et al., 2018).

Spatial distribution data were used to develop a predictive habitat suitability model for the selected kinfisher species. This model was then overlaid with forest function data from Java Island. The 2023 forest function designation map from the Ministry of Environment and Forestry (KLHK) of the Republic of Indonesia was used in the analysis. Indonesia classifies its land into forest areas and area for other land uses (APL). Predicted habitat suitability values at field occurrence locations for kingfisher species were extracted from the habitat suitability maps using the extract values to points tool in ArcGIS. These values were summarised using boxplots to evaluate the correspondence between model predictions and species occurrence records.

Prior to model development, multicollinearity diagnostics were conducted on ten independent environmental variables. The dendrogram of the Pearson correlation analysis indicated high multicollinearity for several variables, including elevation, temperature, radar VH and radar VV (Figure 3). Specifically, elevation and temperature suggested strong interdependence. This result is consistent with established ecological theory, in which temperature gradients are closely linked to elevation changes, particularly in mountainous regions (Loewen et al., 2023). Given this relationship, and to minimize collinearity, elevation was retained in the model as a proxy variable.

Additionally, PALSAR HV, PALSAR HH and NDVI (Normalized Difference Vegetation Index) were found to be perfectly collinear. To address this issue, only NDVI was retained for further analysis, given its strong ecological relevance as a proxy for vegetation density. Consequently, six independent variables were selected for inclusion in the random forest (RF) model and further evaluated using a variable importance analysis (varImp). The variables were built-up land (built_up), elevation (elevation), vegetation density (NDVI), distance to the nearest water body (dist_water), rainfall (precipitation), and vegetation structural complexity (radar VH). All models were developed using the six independent predictor variables. The RF analysis resulted in differing proportions of predicted suitable habitats across the nine species.

The model achieved AUC-ROC values exceeding 0.70 for all nine species, indicating strong predictive capacity (AUC range = 0.75–0.93; mean = 0.83; SD = 0.06). Repeated 5-fold cross-validation results were consistent across species, ranging from 0.74 ± 0.01 to 0.93 ± 0.01, suggesting stable model performance. Sensitivity values ranged from 0.68 to 0.86, while specificity ranged from 0.67 to 0.86. The optimal classifications thresholds derived from the ROC analysis varied slightly among species, ranging from 0.47 to 0.63. These findings demonstrate that the habitat suitability models performed well overall, although model accuracy varied depending on the level of presence data available for each species.

Based on the varImp values, three key variables were identified, with elevation and precipitation emerging as the main factors influencing the distribution of most kingfisher species in Java (Table 2). However, for the white-throated kingfisher, rufous-backed kingfisher, Javan blue-banded kingfisher and blue-eared kingfisher,the variable dist_water influenced the presence of suitable habitat.

To explore the functional relationships between the independent variables and species presence probabilities, the itsdm package in R was employed to generate response curves (Figure 5). These visualizations revealed species-specific variation in how predictor variables influenced the likelihood of occurrence. For example, increasing elevation was generally associated with a decline in the probability of occurrence for several species, namely the cerulean kingfisher, sacred kingfisher, white-throated kingfisher, and banded kingfisher. The probability of the presence of thecerulean kingfisher and sacred kingfisher declinedwhen precipitation exceeded > 1750 mm.

In contrast, elevation showed a weak and non-linear relationship with the presence of the Javan blue-banded kingfisher, as indicated by a relatively flat response curve centered around a probability of ~0.4. This suggests either limited variation in the elevation data for this species or the existence of complex interactions with other environmental variables. For the Javan blue-banded kingfisher, the probability of occurrence increased notably with proximity to water bodies (within 0–300 meters) and with higher precipitation levels (ranging from 2,575 to 3,250 mm), highlighting the species' preference for moist, riparian habitats. Increased precipitation (> 1750 mm) also increased the presence probability of the blue-eared kingfisherandcollared kingfisher.

The distribution of suitable habitat varies based on forest area function and kingfisher species on Java. For all species, the majority of these suitable habitats were located within areas designated for other land uses (Figure 6). This indicates that, in addition to forest status, environmental variables contribute to habitat suitability. These areas are not only important habitats for resident species, but also for migratory ones, including the sacred kingfisher, which is widely distributed across the island. As shown in (Figure 4), distribution is concentrated along the island's periphery and river channels.

The collared kingfisherand Javan kingfisher demonstrate the most extensive spatial distributions among resident species, reflecting their ecological adaptability and capacity to occupy a wide range of land-use types. In contrast, the Javan blue-banded kingfisher exhibits the smallest predicted suitable habitat area, totaling 998,128.4 hectares (Table 2), but shows the highest proportion of suitable habitats in forest areas (48.78%). The distribution of the Javan blue-banded kingfisher is fragmented, occurring primarily in isolated patches associated with mountainous and valley regions that include more extensive natural landscapes or designated protected forest areas, in contrast to the broader and more anthropogenically influenced distributions observed for other kingfisher species (Figure 5).

The habitat suitability model was validated based on habitat suitability values ​​extracted from occurrence locations (indicated by citizen science data) of each kingfisher species on the island of Java. The boxplot shows that the habitat suitability values ​​for the occurrence for all species have a median value > 0.80 (Figure 7). This indicates that the majority of occurrence records recorded by citizens are located within areas predicted by the model as highly suitable habitats. Although there were several outlier occurrence locations with low suitability values, these ​​represent only a small proportion of the data, where environmental conditions are not fully represented by the dominance of predictor variables.

The distribution of habitat suitability values ​​for all species varies in terms of habitat use (Figure 6) and ecological requirements (Figure 5). In the boxplot, the Javan kingfisher species shows the greatest variability, with outliers ranging from approximatelty 0.33 to 0.50 and with the lowest median suitability value (~0.82), indicating a broader habitat tolerance than the other species. Conversely, the sacred kingfisher, cerulean kingfisher and rufous-backed kingfisher species have high median habitat suitability values (> 0.90), indicating strong agreement in model predictions and occurrences. The Javan blue-banded kingfisher, the species with the smallest habitat suitability area yet the highest percentage within forested areas, exhibits a high median habitat suitability value (~0.92) with relatively small variation. Similar variation in habitat suitability is also observed in the white-throated kingfisher, with outlier values that are not particularly low (~0.75). The overall variation in the model effectively captures the environmental conditions associated with species occurrences in the field.

The ongoing global biodiversity crisis requires rapid assessment of knowledge gaps and the setting of conservation priorities based on large-scale biodiversity data. While earlier studies on kingfisher relied heavily on direct field observations and telemetry methods (Kesler, 2012), our approach reveals broader spatial patterns using accessible data sources. In the study, we modeled the potential distribution of nine kingfisher species across Java Island and examined the relationship between their occurrence and key environmental variables. The findings demonstrate that species distribution models using varied citizen science platform data—analysed with the RF algorithm can provide reliable performance; avoid overfitting in spatial prediction (Valavi et al., 2022); and prevent over-confidence in predictive abilities (Akselrud, 2024). Consequently, model performance was contingent upon repeated cross-validation and several evaluation metrics. The combination of machine learning techniques and massive datasets in citizen science projects is a useful approach to identifying spatial distribution (Martín et al., 2021). The consistency between model outputs and established ecological knowledge for each species reflects species-environment relationships and validates the robustness of the modelling approach (Mahmoud et al., 2024). This allows the model to uncover patterns that reflect the ecological responses of species to environmental gradients. The potential distribution map produced the by RF model therefore provided more accurate and informative results (Zhao et al., 2022). Although natural colonisation of a species takes considerable time, the target species may not always be present or detected (Frakes & Knight, 2021).

The predictive models showed relatively good accuracy for all species, with distinct associations with specific habitat features. Notably, despite Java being the most densely populated island in Indonesia, with extensive non-forest areas, built-up areas did not significantly influence the distribution of any kingfisher species (Table 2). Instead, most species were associated with human-impacted, non-forest habitats (Figure 6), underscoring their generalist tendencies (Ariyanto et al., 2024). However, environmental variables such as elevation and precipitation emerged as primary factors shaping habitat suitability for nearly all species. Furthermore, proximity to water bodies (dist_water) was particularly important for the white-throated kingfisher, rufous-backed kingfisher, Javan blue-banded kingfisher, and blue-eared kingfisher, suggesting species-specific habitat preferences within the group. Citizen science observers often recorded kingfisher in human-modified landscapes, such as rice fields in Java, which are commonly located near water sources. Around 75% of Indonesia’s irrigated rice fields use surface or groundwater, while the rest rely on natural sources such as rain, tidal rivers, seawater or seepage (Tirtalistyani et al., 2022). Kingfisher species are also commonly found in urban areas that maintain or create wetlands and water bodies (Xie et al., 2022). These water features play an important role in supporting bird presence, even in degraded habitats (Paga et al., 2022). Hence, we cannot exclude the possibility that some models may not have performed optimally because unsampled combinations of environmental predictors were not taken into account (e.g., interspecific interactions, key specific resources or local disturbance). Future research which compares their distribution across different landscapes, such as watersheds with varying human population densities, could provide deeper insights into their home range and habitat preferences on a finer scale. Caution in interpreting habitat suitability on the local scale is necessary to ensure the results remain valid in representing the actual occurrence or absence.

Rainfall plays a critical role in shaping the distribution patterns of many kingfisher species on Java, as indicated by the predictive habitat models. Being a large and topographically varied island, Java experiences greater rainfall variability than smaller islands, which contributes to diverse microhabitats across the landscape (De Ruyck & Koper, 2024). However, ongoing climate change is expected to significantly alter rainfall regimes and temperature patterns, which may disrupt the ecological conditions that currently support kingfisher populations (Liang et al., 2021; McPherson et al., 2025). These changes could lead to shifts in species distributions, making some areas that are presently suitable become less hospitable in the future. This is particularly concerning given that the majority of suitable habitats identified for kingfisher species are located in non-forest, human-modified landscapes, which are also under pressure from land-use change and development. As such, understanding the interplay between climate variability—especially rainfall—and species ecology is essential for long-term conservation planning. Future research should integrate climate projections and rainfall variability to improve predictive models and assess potential habitat shifts, thereby helping to anticipate the impacts of climate change on kingfisher species across the Island.

The presence probability of endemic, migratory and protected kingfisher species in Java is strongly influenced by vegetation density, as indicated by the NDVI. For example, the banded kingfisher is typically associated with protected hill forests with dense vegetation and intersected by rivers, highlighting its preference for intact forest ecosystems (Hazwani et al., 2023). Vegetation near water bodies is ecologically important as it provides perching sites used by kingfishers to detect and capture prey (Naher & Sarker, 2014). Similarly, species such as the cerulean kingfisher and sacred kingfishershow higher distribution probabilities along the northern coastal regions of Java (Figure 4), where large mangrove areas exist. However, these habitats are under increasing threat due to human activities such as fish pond development and aquaculture expansion (Hanifa et al., 2024).

Despite variations in vegetation cover, many kingfisher e species continue to associate closely with wetland environments. For instance, the collared kingfisheris commonly found near rivers, even in areas with abundance of vegetation (Biswas & Deuti, 2025), while the cerulean kingfisher inhabits open coastal zones with sparse vegetation (Van Balen, 1998). The lowland riverbanks of the island are often exposed and heavily influenced by human disturbance (Prihestiwi et al., 2023), yet remain vital for feeding and breeding. Rivers provide not only food sources, such as fish, for the collared kingfisher (Biswas & Deuti, 2025), but also nesting sites, including steep soil banks used by theJavan kingfisherfor burrow construction(Taufiqurrahman et al., 2020). Additionally, the blue-eared kingfisheris known to stay along riverbanks during the non-breeding season, suggesting riparian dependence (Palkar, 2016). These findings emphasise that vegetation density—especially in riparian zones—has a stronger influence on kingfisher habitat suitability in Java than vegetation structure or forest cover alone, underlining the importance of conserving water-associated landscapes even in modified environments. We utilised the Ministry of Environment and Forestry's forest area function classification in this analysis by including water bodies covering 0.001% of the total area of ​​Java Island. The modeling results indicate that water bodies constitute

a very small percentage (0 - 0.01 %) of the suitable habitat area for each kingfisher species. Given the very small proportion of water bodies at the geographic scale on the island, they were retained in the analysis but combined with APL to facilitate categorization based on forest area function. Kingfisher species have strong ecological associations with water areas and tree cover (Barik et al., 2022; Shifa et al., 2023). Therefore, further research on finer scales is needed to understand the mechanisms driving kingfisher habitat selection at the microsite scale, particularly with regard to water body characteristics and riparian vegetation. Riparian ecosystems are widely surrounded by landscapes that have been modified by human activities such as intensive agriculture and artificial water infrastructure to support agriculture (Petit et al., 2023), as is the case on Java.

The predicted habitat distribution of kingfisher species provides important information to support conservation planning at local scales by identifying suitable habitats based on key environmental variables in managed and unmanaged areas. Many of these habitats are located outside designated conservation zones, underscoring the need to expand protection measures to non-conservation areas across Java. This is particularly critical for rare species such as the Javan blue-banded kingfisher, which require connectivity between suitable and fragmented habitats to maintain viable populations. Therefore, preventing land-use conversion and establishing ecological corridors to connect habitat patches should be prioritised to facilitate species movement and reduce isolation (Ariyanto et al., 2024; Imron et al., 2023). Although there is limited evidence of hunting pressure on kingfisher species, threats from land conversion and declining environmental quality—particularly water pollution—pose more significant risks to their survival (Barik et al., 2022; Shifa et al., 2023). Further investigation of the potential roles of markets in their distribution might provide clearer conservation intervention (Imron et al., 2025). Active community involvement through citizen science initiatives on the island of Java also plays a crucial role not only in expanding species occurrence data (Winnasis et al., 2018), but also in enabling long-term monitoring of population trends and distributions, especially for endemic, migratory, and protected species (Squires et al., 2021).