An Integrated Spatio-Temporal Analytics Framework for Smart Geotourism: Combining K-Means Clustering and SARIMAX-Based Visitor Demand Forecasting in Bayah Dome Geopark

Modern geotourism management faces the dual challenge of anticipating temporal fluctuations in visitor demand and managing a spatially diverse portfolio of attractions. Bayah Dome Geopark illustrates this challenge through its 49 geosites, which differ in accessibility, infrastructure, and experiential characteristics. Integrating visitor demand forecasting with geosite segmentation provides a coherent framework for scenario-based destination planning [1, 2].

Previous studies have generally addressed tourism demand forecasting and destination segmentation as separate analytical tasks [3, 4]. Forecasting commonly focuses on aggregate demand, whereas clustering provides static spatial profiles without accounting for time-series fluctuations. This separation limits the ability to interpret forecasted demand alongside specific geosite characteristics. In tourism demand forecasting, the deliberate selection of exogenous variables and the integration of contemporaneous data represent pivotal factors for building operational forecasting models [5]. To support tourism planning and destination development, stakeholders can adopt analytical techniques such as cluster analysis to identify meaningful groups of destinations based on shared characteristics [6]. This persistent separation between temporal and spatial analysis represents a significant research gap, limiting the development of proactive and coordinated strategies for managing visitor flows in real-world operational contexts.

This study addresses this gap by integrating K-Means clustering and Seasonal Autoregressive Integrated Moving Average with Exogenous Regressors (SARIMAX) forecasting within a scenario-based decision-support prototype. K-Means identifies interpretable profiles across 49 geosites, while SARIMAX forecasts geopark-level daily visitor demand using weekly seasonality and external variables. The integration illustrates how temporal forecasts and spatial profiles could conceptually support management and trip-planning scenarios without claiming demonstrated behavioral or operational impacts.

This study adopts a structured Research and Development (R&D) methodology. This approach was selected to reflect the study's dual objective: generating novel analytical insights and translating these findings into a tangible, operational smart geotourism application. The R&D process comprised the following phases:

Preliminary research: This initial phase involved a comprehensive literature review to identify the research gap, followed by the collection of time-series and cross-sectional data, and extensive data preprocessing to ensure data quality and suitability for analysis.

Model development: The analytical core of the research, this stage focused on constructing the two complementary machine learning pipelines: spatial segmentation using K-Means clustering and temporal forecasting using a SARIMAX model.

Product development: This stage involved translating the outputs from the analytical models into a functional, web-based application with distinct, user-centric interfaces for both geopark management and end-user tourists.

Validation: The final phase consisted of a dual evaluation. First, the statistical accuracy of the predictive models was assessed against established metrics. Second, the technical feasibility and decision-support potential of the integrated framework were evaluated through a scenario-based proof-of-concept application. The latter evaluation did not represent real-world user validation or field deployment.

The first pillar, geosite segmentation, addresses the types and spatial distribution of attractions exist and 'where' they are concentrated by partitioning the 49 geosites into interpretable groups. The second pillar, visitor forecasting, tackles the operational question of 'when' visitor demand will fluctuate by modeling daily arrivals and incorporating exogenous variables. The true novelty lies not in the isolated application of these techniques, but in their deliberate synergy throughout the R&D lifecycle. This integrated foundation supports a transition toward proactive, data-driven geotourism management. The entire methodological pipeline is visually summarized in Figure 1.

tu_pian_40.png

Figure 1. Research framework

3.1 Data and pre-processing

The analysis was based on two distinct datasets curated to support the core analytical pillars of this research. The first was a daily time series dataset for visitor forecasting, spanning from January 2023 to December 2024 (731 observations). Daily visitor counts across 49 individual geosites were aggregated by date to construct a unified, geopark-level time series. This target variable was supplemented by exogenous regressors selected for their potential influence on geotourism demand, including holiday indicators, online search interest, and morning weather conditions [20]. The second was a cross-sectional dataset detailing the unique attributes of the 49 geosites for the segmentation task. These features represented key dimensions of tourist experience and operational management, including geosite type, accessibility, visitor appeal, and level of development. The variables used to characterize the geosites for the K-Means clustering process are summarized in Table 1.

Table 2 presents the variables included in the visitor forecasting model, comprising the dependent variable and exogenous regressors.

Table 1. Variables used in clustering with types and description

Table 2. Variables used in forecasting, with types and descriptions

The exogenous variables were aligned by date with the daily geopark-level visitor series and comprised a binary Holiday Indicator, Google Trends Index, mean morning cloud cover, and mean morning relative humidity (RH). The Holiday Indicator distinguishes “Holiday” and “Non-Holiday” dates, while the Google Trends Index represents normalized relative search interest. Meteorological variables represent the corresponding mean morning measurements. For the hold-out evaluation, these exogenous observations were aligned with the respective test dates. Because this study does not forecast future Google Trends or weather values, the reported SARIMAX results represent conditional forecasting performance based on the prepared exogenous inputs rather than a fully autonomous forecasting system.

3.2 K-Means clustering

To perform geosite segmentation, this study employed K-Means clustering to identify homogeneous groups of geosites based on their shared characteristics [21, 22]. The selection of clustering features was guided by domain relevance and prior literature on geotourism destination evaluation. Specifically, features were chosen to represent three key dimensions critical for visitor flow management: accessibility (e.g., travel time, road condition), attractiveness (e.g., visitor ratings, review counts, thematic attributes), and operational capacity (e.g., facilities, entrance fee). This feature selection strategy was designed to ensure that the resulting clusters were interpretable and directly relevant for destination management [23, 24].

The final feature set entered into the K-Means analysis consisted of geosite type, distance, travelling time, visitor rating, road condition, entrance fee, number of facilities, photography, education, adventure, family recreation, and rock type. The ‘Geosite Name’ attribute was excluded from the clustering feature matrix and retained solely as an identifier.

Prior to clustering, all numerical features were standardized using Z-score normalization to prevent variables with larger numeric ranges from disproportionately influencing the distance-based K-Means algorithm. Categorical variables were encoded using dummy variables to allow their inclusion in the clustering process. The optimal number of clusters was determined using the Elbow Method, which evaluates the trade-off between cluster compactness and model simplicity. The resulting within-cluster variation across candidate values of k is presented in Figure 2.

A clear inflection point was observed at k = 4, indicating that four clusters provide a parsimonious and interpretable segmentation of the geosite dataset.

tu_pian_41.png

Figure 2. Elbow method for selecting k in K‑Means

To further assess the clustering solution, internal clustering validity and stability across random initializations were evaluated. For the selected four-cluster solution, the silhouette coefficient was 0.1751, the Davies–Bouldin index was 1.6113, and the Calinski–Harabasz index was 10.314. These metrics yielded complementary insights into cluster cohesion, separation, and compactness. Cluster stability was evaluated across 50 random K-Means initializations using the Adjusted Rand Index (ARI). The mean pairwise ARI across the 50 runs was 0.631 (SD = 0.153), with values ranging from 0.369 to 1.000, while the mean ARI relative to the final solution using seed = 42 was 0.671. These results demonstrated a moderate level of consistency in cluster assignments across varying initializations.

The resulting four clusters were subsequently interpreted based on their mean feature profiles. Cluster 0 was characterized as Accessible Cultural & Bio-Diversity Geosites, Cluster 1 as Nearby but Underdeveloped Geological Geosites, Cluster 2 as Developed Geo-Adventure & Family Hubs, and Cluster 3 as Distant but Highly-Rated Geo-Adventure Sites. These profiles summarize differences in accessibility, infrastructure, geosite type, and visitor-oriented characteristics across the four segments.

As a descriptive comparison, one-way Analysis of Variance (ANOVA) was conducted on the numerical features included in the clustering process, namely distance, travelling time, visitor rating, entrance fee, and number of facilities. Since these variables were directly used to construct the K-Means solution, the resulting between-cluster differences were not treated as independent evidence of cluster validity. ANOVA assumptions were verified, and the Kruskal–Wallis test was additionally used as a non-parametric robustness check where appropriate. The results showed significant differences across clusters for distance, travel time, entrance fee, and facility count, whereas visitor rating exhibited no significant variation among the segments. These descriptive comparisons help characterize the key attributes distinguishing the four segments, while the internal validity and stability metrics provide the primary evidence for the quality of the clustering solution.

The identified cluster profiles provide a strategic foundation for differentiated destination management, including infrastructure prioritization, targeted promotion, and visitor distribution planning across the Bayah Dome Geopark. Within the prototype web application, these profiles are integrated to generate tailored recommendation scenarios aligned with specific visitor preferences such as adventure, photography, education, and family recreation. However, these applications represent scenario-based decision support and have not yet been validated through field deployment or real-world user testing.

3.3 SARIMAX forecasting

To forecast daily visitor arrivals, a SARIMAX model was employed owing to its capability to simultaneously capture historical demand dynamics, weekly seasonality, and exogenous explanatory variables [25]. The model was developed using the daily geopark-level visitor time series introduced in Section 3.1. The forecasting procedure followed a structured pipeline involving chronological data splitting, stationarity assessment, model identification, parameter estimation, and out-of-sample evaluation.

To evaluate forecasting performance under realistic operational conditions, the dataset was chronologically divided into training (80%, n = 584) and testing (20%, n = 147) subsets. The training dataset was used for parameter estimation and model selection, whereas the testing dataset was reserved exclusively for out-of-sample forecasting evaluation. A chronological split was adopted to preserve temporal dependencies and prevent information leakage between training and testing data. First, the stationarity of the daily visitor time series was examined using the Augmented Dickey–Fuller (ADF) test [26]. The ADF test indicated stationarity at the 5% significance level (p < 0.05). Nevertheless, the ADF result was not used as the sole criterion for establishing the differencing orders. Instead, regular and seasonal differencing orders were treated as model-selection parameters and systematically evaluated within the SARIMAX grid search framework. This approach enabled candidate models with differing combinations of non-seasonal and seasonal differencing to be compared via information criteria while accounting for the observed weekly seasonality. Seasonal patterns were further evaluated through time-series decomposition and autocorrelation diagnostics.

Model order identification was initially diagnosed through visual inspection of the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots. To systematically evaluate alternative differencing and autoregressive-moving-average structures, a grid search involving 120 candidate configurations was conducted across non-seasonal and seasonal parameters, given a weekly seasonal cycle (s = 7). Candidate models were evaluated using the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) for model selection. The selected SARIMAX specification, SARIMAX (0,1,2)(0,1,2,7), achieved the lowest AIC (3727.267) and BIC (3766.202).

The selected differencing orders should therefore be interpreted as the outcome of a comparative model selection process rather than as a direct implication of the ADF test. In particular, the final specification was retained because the combination of regular differencing (d = 1), seasonal differencing (D = 1), and the corresponding MA terms provided the best information-criterion performance among all evaluated candidate configurations.

In addition to historical visitor counts, exogenous variables were incorporated to represent contextual factors potentially associated with visitor demand. These variables consisted of holiday indicators, Google Trends search interest, morning cloud cover, and RH. Although the estimated coefficients exhibited varying levels of statistical significance, these variables were retained in the final specification based on their theoretical relevance and contextual explanatory value. Specifically, calendar events, meteorological factors, and online search behavior provide valuable signals for capturing visitor-demand fluctuations in an outdoor geotourism context.

To assess the incremental forecasting value of the exogenous variables, an additional benchmark experiment was conducted by comparing the selected SARIMAX model with an equivalent SARIMA model without exogenous regressors. Both models utilized the same autoregressive, differencing, moving-average, and seasonal orders, alongside identical chronological training and testing subsets, to ensure a consistent comparison. Model fit was evaluated using AIC and BIC, while out-of-sample forecasting performance was assessed via Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and Symmetric Mean Absolute Percentage Error (sMAPE).

Model adequacy was further assessed through residual diagnostic tests, including the Ljung–Box test for residual autocorrelation, the Jarque–Bera test for residual normality, and the Autoregressive Conditional Heteroskedasticity Lagrange Multiplier (ARCH-LM) test for conditional heteroskedasticity. The complete diagnostic test outputs are detailed in the Results section, while the estimated coefficients, standard errors, z-statistics, and p-values of the selected SARIMAX model are presented in Table 3.

Table 3 summarizes the estimated parameters of the selected SARIMAX specification. The MA(1) and seasonal MA(7) components were statistically significant, whereas the exogenous variables did not reach individual statistical significance. Nevertheless, the exogenous variables were retained as theoretically motivated contextual predictors and subsequently benchmarked against an equivalent SARIMA baseline model without exogenous inputs. The seasonal MA(7) coefficient was estimated at −0.9986, which lies close to the conventional invertibility boundary. Although the parameter was statistically significant, its near-boundary estimate warrants cautious interpretation and may indicate that the selected specification is sensitive to the strong weekly seasonal structure in the visitor series. This potential numerical limitation was considered when interpreting the residual diagnostics and forecasting performance. The temporal dependence and seasonal structure of the visitor series were further examined using the ACF and PACF plots presented in Figure 3.

Table 3. Estimated parameters of the final SARIMAX model

tu_pian_42.jpg

Figure 3. Autocorrelation (ACF) and Partial Autocorrelation (PACF) plots

Figure 3 presents the ACF and PACF diagnostics used to examine the temporal structure of the visitor data. In the ACF plot (top), significant and recurring correlation spikes are observed at lags that are multiples of 7 (lags 7, 14, 21, and so on), indicating a weekly seasonal pattern (s = 7). Meanwhile, the PACF plot (bottom) provides preliminary information for determining appropriate autoregressive and moving-average orders. The shaded region represents the 95% confidence interval, with observations outside this region indicating statistically significant autocorrelation.

3.4 Integration logic

The integration logic connects the temporal and spatial analytical components of the framework. SARIMAX provides geopark-level visitor-demand forecasts, while K-Means generates geosite profiles based on accessibility, infrastructure, and visitor-oriented characteristics. Together, these outputs could support scenario-based consideration for promotion, resource allocation, and alternative geosite planning during anticipated periods of high demand. The framework is intended to provide decision-support information rather than evidence of actual visitor redistribution or behavioral change.

3.5 Proof-of-concept validation of the integrated framework

The integrated framework was evaluated through a scenario-based proof-of-concept to assess its technical feasibility and intended decision-support logic rather than actual changes in visitor behavior. A representative high-demand holiday scenario was constructed using the SARIMAX forecast, while K-Means profiles were deployed to identify geosites with different accessibility, facility, and visitor-oriented characteristics for potential alternative planning.

The scenario illustrates how forecast information and geosite profiles can support data-driven decision-making in strategic promotion, resource allocation, and visitor trip planning. However, it does not provide empirical evidence of changes in tourist decisions, congestion, or visitor redistribution. Such effects require future field deployment, usability testing, behavioral observation, and longitudinal evaluation.

This study developed a scenario-based framework integrating K-Means clustering and SARIMAX forecasting for spatial geosite segmentation and temporal visitor-demand analysis. K-Means produced four interpretable geosite profiles, while SARIMAX achieved an sMAPE of 6.34% on the hold-out test set. The 95% prediction intervals achieved 97.28% empirical coverage; however, mild autocorrelation at lag 7, substantial residual non-normality, and systematic underestimation of holiday peaks indicate that the forecasting outputs should be interpreted as useful predictive performance rather than full statistical validation.

The integrated outputs were deployed within a proof-of-concept web application designed to support promotional planning, resource allocation, and trip-planning scenarios. However, the prototype does not provide empirical evidence of changes in tourist behavior, site congestion, visitor redistribution, or management effectiveness.

Future research should evaluate the framework through field deployment, usability testing, and longitudinal observation. Additionally, subsequent studies could explore advanced forecasting models, alternative clustering techniques, and additional real-time data sources.

All authors contributed to the preparation and completion of this study, including research planning, literature review, system design, system testing, data analysis, manuscript preparation, and final revision.