Communications in Science and Technology 10. 135Ae147 COMMUNICATIONS IN SCIENCE AND TECHNOLOGY Homepage: cst. Leveraging machine learning and open accessed remote sensing data for precise rainfall forecasting Bambang Kun Cahyonoa,*. Muhammad Hidayatul Ummaha. Ruli Andarua. Neil Andikab. Adjie Pamungkasc. Hepi Hapsari Handayanid. Paramita Atmodiwirjoe . Rory Nathanf Geodetic Engineering Department. Universitas Gadjah Mada. Yogyakarta 55281. Indonesia Civil and Environmental Engineering Department. Universitas Gadjah Mada. Yogyakarta 55281. Indonesia Urban and Regional Planning Department. ITS. Surabaya 60117. Indonesia Geomatics Engineering Department. ITS. Surabaya 60117. Indonesia Department of Architecture. Universitas Indonesia. Depok 16425. Indonesia Department of Infrastructure Engineering. University of Melbourne. Parkville 3010. Australia Article history: Received: 11 January 2025 / Received in revised form: 10 June 2025 / Accepted: 13 June 2025 Abstract Rainfall forecasts are essential for human activities enabling communities to anticipate any impacts. Rainfall events correlate with other natural and hydro-meteorological phenomena, which can be used in modeling and prediction. This study used daily CHIRPS for the Gajahwong watershed in Yogyakarta. Indonesia as the precipitation data. It also used Sea Surface Temperature. Land Surface Temperature (Day and Nigh. Minimum and Maximum Temperatures. Solar Radiation. Wind Speed (U and V component. Cloud Pressure (Top and Bas. , and Cloud Height (Top and Bas. as the parameters. Further, data processing was performed by means of the Google Earth Engine (GEE) platform. Machine learning methods, including Support Vector Regression. Gradient Boosting Regression. Random Forest, and Deep Neural Networks, were The correlation analysis revealed that only the Wind Speed V-component showed significant correlation with rainfall, other seven parameters showed moderate and four showed weak ones. Meanwhile, accuracy assessments indicated that Support Vector Regression had the most accurate predictions accompanied by Root Mean Square Error (RMSE). Mean Absolute Error (MAE). Mean Squared Error (MSE). R2, and Coefficient Correlation (CC) at 1. 366, 0. 947, 1. 866, 0. 948 and 0. 982 respectively. This study demonstrated that utilizing openly accessible atmospheric datasets processed through the GEE could yield reliable rainfall predictions, facilitating informed decisions on a wide scale. The methodology is adaptable and can be reproduced for any comparable research or operational purposes. Keywords: Precipitation prediction. rainfall forecast. machine learning. Gajahwong River. hydro-meteorological big data Introduction Rainfall is a primary factor of the hydrological cycle, essential for sustaining ecosystems and maintaining environmental balance on Earth . The process begins with evaporation, followed by condensation, cloud formation, and precipitation . Water, once on the ground, undergoes interception, infiltration, transpiration, and evaporation, completing a hydrological cycle . Rainfall distribution in a specific area is determined by various atmospheric factors, including surface temperature, humidity, air pressure, wind speed, cloud cover, and solar radiation . In such condition, satellite imagery serves as a reliable tool for detecting and documenting atmospheric phenomena and conditions . As a tropical country. Indonesia * Corresponding author. Tel. : 62-274-520226 . fax: 62-274-520226 Email: bambangkun@ugm. https://doi. org/10. 21924/cst. experiences substantial rainfall with some equatorial regions receiving it year-round, while others have distinct seasonal patterns . Weather conditions significantly impact human activities such as in agriculture . planting, seedling, maintenance, and harvestin. , transportation . and and ai. , tourism . utdoor and natural attraction. , flood mitigation . , and water resource management . For this, accurate short-term and long-term weather forecasts become vital across multiple sectors . Accurate weather forecasting is deemed vital for anticipating any events that can disrupt daily activities . Climate change has further altered rainfall patterns, making precise prediction increasingly important . Beyond mere delays in human activities, extreme weather events such as prolonged droughts, flash floods, infrastructure damage, and landslides pose serious threats to many regions . To mitigate these negative impacts, reliable rainfall This open access article is distributed under a Creative Commons Attribution (CC-BY) 4. 0 license Cahyono et al. / Communications in Science and Technology 10. 135Ae147 prediction and modeling are crucial. Advanced forecasting models provide early warnings . , helping to minimize any risks to life, property, infrastructure, and land . While communities have historically relied on local knowledge and wisdom to predict rainfall . , the growing impact of climate change on weather conditions . necessitates more precise and technologically advanced approaches. Accurate rainfall predictions benefit various sectors by enabling efficient water resource management and effective disaster mitigation. Leveraging big data from diverse weather parameters and applying machine learning algorithms can significantly enhance prediction models . ,17,. Techniques such as classification and regression allow for the analysis of historical data to forecast future rainfall patterns . Various methods for rainfall prediction have been developed, including both statistical and machine learningbased approaches . Statistical modeling derives equation models based on existing data . ata-drive. , utilizing techniques such as Simple Regression Analysis (SRA). Decomposition. Exponential Smoothing (ES). Autoregressive Integrated Moving Average (ARIMA) . , and Least Squares Adjustment (LA) . Meanwhile, machine learning-based rainfall models such as Neural Networks (NN). Random Forest (RF). Gradient Boosting (GB). Support Vector Machines (SVM) . K-Nearest Neighbors (KNN), and Genetic Programming (GP) . are considered highly accurate in identifying patterns in existing data. Machine learning algorithms are able to effectively identify rainfall patterns from historical event data, enabling prediction and regression. These approaches are categorized as univariate forecasting, which relies on patterns in historical data involving a single variable . Univariate forecasting is advantageous in consideration to its simplicity in interpretation and computational efficiency . In contrast, multivariate time series forecasting incorporates multiple variables as the inputs for predicting future rainfall . This approach is believed to enhance accuracy, capture complex interrelationships among variables, address external factors, and improve robustness . Machine learning in this context analyzes correlations between various hydrological and hydro-meteorological phenomena and rainfall occurrences. These correlations can then be used as parameters that influence modeling and prediction . Atmospheric hydro-meteorological conditions correlated with rainfall include evaporation, sunshine, wind speed, humidity, cloud properties, temperature . , dew point, wind direction, visibility . Southern Oscillation Index. NINO 3. 4 Index . Madden-Julian Oscillation (MJO). Northern Oscillation Index (NOI), and Quasi-Biennial Oscillation (QBO) . Hydro-meteorological data used as covariates are recorded by various sensors attached on satellite platforms. This is an effective alternative to address the limitations of ground station observation data in terms of quantity and accessibility . The quality of satellite data, however, is significantly determined by the resolution of the recording sensor . , including spatial resolution . etail of coverag. , temporal resolution . requency of observation. , and spectral resolution . ensor capabilitie. With the decades of recorded observations, historical events can be analyzed to uncover patterns and trends . ,25,. This wealth of satellite data can be considered "big remote sensing data," ready for processing and analysis to address various challenges . Relevant institutions provide a wide range of hydroclimatic data freely accessible to public. These datasets are collected through continuous observations by satellites, ground stations, or a combination of both . Moreover, these data can be accessed and processed through free platforms such as Google Earth Engine (GEE) . ,30,. Many researchers are widely using GEE, which offers multivariate data with strong correlations to rainfall patterns . Based on the current conditions and existing literature, the integrated use of complex atmospheric variables and multivariate approaches in rainfall modeling, so far, remains significantly underexplored. Similarly, the application of the Google Earth Engine platform for both univariate and multivariate rainfall prediction using machine learning techniques seems still limited. These gaps highlight a valuable opportunity to investigate the importance of accurate rainfall forecasting and the potential of leveraging publicly available hydro-meteorological data. This study, in turn, aims to forecast rainfall over the next five epochs using machine learning-based univariate and multivariate models. The novelty of this study lies in the use of complex atmospheric parameters such as Sea Surface Temperature. Land Surface Temperature (Day and Nigh. Minimum and Maximum Temperatures. Solar Radiation. Wind Speed (U and V component. Cloud Pressure (Top and Bas. , and Cloud Height (Top and Bas. in rainfall modeling and prediction via the GEE platform. Each parameter underwent a correlation analysis to determine its influence on rainfall and alignment with rainfall patterns. Additionally, multiple machine learning methods were evaluated and compared to identify which approach delivered the most accurate predictions. Materials and Methods Area study This research was conducted in the Gajahwong River watershed, which spans from the peak of Mount Merapi to Bantul Regency in the Special Region of Yogyakarta. Indonesia. The area is located between coordinates . A32'26. 72"S, 110A26'45. 52"E) at the upper part of Mount Merapi and . A50'21. 41"S, 110A23'47. 35"E) in the downstream The watershed features a steep river channel with significant topographic changes in the northern area, transitioning to low-lying, nearly flat land in the southern area that gradually slopes downward. The low-lying areas have a slope of 10%-15% and are the most vulnerable to flooding . The elevation ranges from 2,905 meters sea level height (SLH) at the peak to 118 meters SLH in the southern part, covering a total area of approximately 44. 40 kmA. Fig. 1 illustrates the study area. This study area was selected with a consideration that the river serves as a vital source of irrigation for rice fields and traverses several key landmarks, including state and private universities, as well as a zoo. Additionally, the area along the river is prone to frequent flooding, particularly during the Cahyono et al. / Communications in Science and Technology 10. 135Ae147 period of heavy and high-intensity rainfall . These flooding events are attributed to the river's proximity to a fault zone . eomorphological facto. , rapid land use changes, and inadequate water management practices. The study area lies within a tropical climate zone, characterized by unpredictable rainfall with high annual variability and an average air temperature ranging from 23AC to 26AC. modeling and prediction, recorded concurrently with the rainfall data. These datasets included sea surface temperature (SST) . , wind speed component-u (WSU), wind speed component-v (WSV), minimum temperature (STMi. , maximum temperature (STMa. , surface net solar radiation (SSR) . , land surface temperature during the day (LSTD), land surface temperature at night (LSTN), cloud top pressure (CTP), cloud top height (CTH), cloud base pressure (CBP), and cloud base height (CBH) . Table 1 depicts the detailed information about the data sources, spatial units, and value units for each dataset. Method Fig. Overview of the Gajahwong River watershed study area Data Preprocessing data Table 1. Detailed information about the dataset ources, spatial resolution, and units of measuremen. Data Precipit SST Data sources CHIRPS https://developers. com/earthengine/datasets/catalog/ UCSBCHG_CHIRPS_DAILY NOAA Optimum Interpolation Sea Surface Temperature (OISST) https://developers. com/earthengine/datasets/catalog/NOAA_CD Spatial Resolution Units 5 KM 0,25 Arc Degree AC WSU WSV R_OISST_V2_1 European Centre for Medium-Range Weather Forecasts (ECMWF) 11 KM 11 KM STMin STMax SSR https://developers. com/earthengine/datasets/catalog/ECMWF_E RA5_DAILY 11 KM 11 KM 11 KM AC J/m^2 LSTD MODIS https://developers. com/earthengine/datasets/catalog/modis 1 KM AC 1 KM 1 KM 1 KM 1 KM LSTN CTP CTH CBP CBH This study employed four machine learning algorithms across basic learning, ensemble learning, and deep learning categories for multivariate rainfall prediction modeling. The selected algorithms included support vector regression (SVR) representing basic learning, random forest (RF) and gradient boosting regressor (GBR) representing ensemble learning, and deep neural network (DNN) representing deep learning. The research was carried out in three key stages: data preprocessing, multivariate rainfall prediction modeling and evaluation, and future rainfall prediction across upcoming epochs. Sentinel 5P https://developers. com/earthengine/datasets/catalog/sentinel-5p This study used daily rainfall data from the Climate Hazards Group Infrared Precipitation with Station (CHIRPS) dataset, spanning the period from 2018 to 2023 . For accuracy, the CHIRPS rainfall data have been calibrated against ground station observations within the study area. Twelve additional datasets with daily temporal resolution were incorporated for In this stage, several preprocessing steps were performed, including feature selection, dataset generation, data scaling, and rainfall data filtering. Feature selection aimed to identify variables significantly correlated with rainfall. The Pearson correlation coefficient, as defined in Eq. , was utilized to measure the relationship between each variable and rainfall. Variables with an absolute correlation value greater than 0. indicating a moderate to strong relationship . , were selected for inclusion in the rainfall prediction process. The parameters description includes yc as the correlation. ycuycn as the value of dataset 1 at point ycn. ycuI as the mean of dataset 1. ycycn as the value of dataset 2 at point ycn. and ycI as the mean of dataset 2. yee= Oc. eoyeOOeyeo I). eoyeOOeyeo I) oc. eoyeOOeyeo I)ya oc. eoyeOOeyeo I )ya . Dataset generation involved defining the model's predictor variables . and target labels . For each variable at a given epoch ycuyc , the predictor variables used in this study included data from the five preceding epochs, denoted as ycuycOe1 , ycuycOe2 , ycuycOe3 , ycuycOe4 , ycaycuycc ycuycOe5 . All variables with a correlation greater than 0. 4 with precipitation, including precipitation itself, were included as the predicting variables. Thus, the total number of predicting variables in this study was ycu y 5, where ycu is the number of variables with a correlation The target label . is defined as the precipitation at a specific epoch, represented as ycyc . Each variable in the dataset had a different range or scale of values, spanning from single units to thousands. To reduce potential biases arising from these differences in value dimensions, this study then applied the Min-Max Scaler normalization technique, rescaling the values of all variables to a uniform range between 0 and 1 . The Min-Max scaling Cahyono et al. / Communications in Science and Technology 10. 135Ae147 process is represented in Eq. , where ycuycnA is the scaled value of the i-th data point, ycuycn is the original value, ycuycoycnycu is the minimum value of the variable, and ycuycoycaycu is the maximum value . yeo Oeyeo yeoAyeO = yeo yeO OeyeayeOyea yeayeCyeo yeayeOyea In addition to preprocessing the prediction variables, this study applied the label data preprocessing. It involved filtering the data to ensure a clean time-series, and reduce the To achieve these purposes, the Butterworth filter, a widely utilized technique in digital signal processing, was The Butterworth filter is known for its maximally flat passband response and is used to remove any undesired frequencies and noise . The filtering process involved data transformation from the time domain to the frequency domain and the application of the Butterworth transfer function. The transfer function is defined in Eq. , where yce represents the frequency at the ycn-th data point, yceyca is the cut-off frequency, and ycA is the filter order . For this study, a low-pass Butterworth filter was applied with a cut-off frequency of 0. 5 and a filter order of 5. The filtering process was implemented using the NumPy library, specifically utilizing the butter function. ya ycyenyeiyeiyeIyeeyeoyeayeeyeiyeO = ya . eN/yeN )yaycA yeE Fig. Filtering results using the Butterworth filter . ata source: CHIRPS accessed and processed using GEE, by the author. Fig. presents the results of rainfall data filtering using the Butterworth filter. Visually, the differences between the original and filtered data were minimal, with only slight deviations observed. The filtering process successfully reduced a number of outliers, resulting in cleaner data better suited for modeling in subsequent steps. Fig. 3 depicts the residual values, representing the differences between the original and filtered data. These residuals ranged from -21. 71 mm/day to 22. 51 mm/day. Most residuals clustered close to zero with an average value of -0. A 3. 62 mm/day. When expressed in absolute terms, the mean residual became 2. 24 A 2. 85 mm/day, indicating that the filtering process effectively preserved the primary data trends while smoothing out noise. The pattern of daily rainfall fluctuations was still observed, but the seasonal pattern was refined by applying Butterworth filtering. Under these conditions, the prediction algorithm was able to learn the timeseries pattern well . Multivariate rainfall prediction modeling and evaluation The prediction dataset was divided using a sequential splitting method with an 80:20 ratio. This meant that 80% of the data were allocated for training, while the remaining 20% was reserved for testing. The training dataset was used to develop the prediction model, and the testing dataset was utilized to evaluate the model's performance. In this paper, the multivariate prediction model was constructed using four machine learning algorithms: SVR. RF. GBR, and DNN. Support Vector Regression (SVR) refers to a regression variant of the Support Vector Machine (SVM) algorithm. SVR was developed by Vapnik . , and is based on statistical and mathematical principles establishing relationships between a set of independent variables and a dependent variable . In SVR, non-linear relationships are addressed by projecting the data into a higher-dimensional feature space where a linear function is used to approximate the This function, represented as a vector, includes an epsilon value to account for uncertainty within the vector SVR employs a deterministic optimization approach to minimize errors. The general formulation of SVR is shown in Eq. , where w is the weight vector in the feature space. I is the transformation function that linearizes the input data in the new feature space, and b is the bias term . SVR offers several advantages, including its effective handling of multidimensional data and relatively low computational requirements . yeo = ycyc Oo yys. yeE, yeo OO ycyeI , yys. OO ycyeI , yeE OO yc Fig. Residual values from rainfall filtering compared to the original data . ata source: CHIRPS accessed and processed using GEE, by the researcher. Random Forest (RF), introduced by Breiman in 2001 . , is an ensemble machine learning model designed for both classification and regression tasks . It operates by constructing multiple decision trees (DT. and utilizing a bagging technique, or bootstrap aggregation, to generate diverse datasets through randomization strategies . Since RF consists of numerous DTs, the final prediction is obtained by averaging the outputs of all individual trees, as described in Eq. Here, yc represents the final output, ycyaycNycn is the output of the i-th decision tree, and ycA is the total number of decision trees generated . The primary strengths of RF include its Cahyono et al. / Communications in Science and Technology 10. 135Ae147 non-parametric nature, high predictive accuracy, robustness against noisy or overfitting data . yeo = ycA OcycA yeO=ya yeoycycyeO Gradient Boosting Regression (GBR), introduced by Friedman in 2001 . , is an ensemble learning model that builds multiple decision trees (DT. iteratively and sequentially . Each subsequent DT is trained to minimize the residual errors of the previous ones. The predicted value of GBR using n decision trees is expressed in Eq. , where yaycNyco represents a weak learner, typically a single decision tree with low individual performance, and yuyco is a scaling factor applied to the tree to minimize residuals. To achieve this. GBR employs gradient descent to adjust the model by updating the predictions based on the residuals of prior estimates . The final model combines the initial estimate with appropriately weighted corrections from subsequent trees. GBR offers notable advantages, including the ability to model complex, non-linear relationships and a relatively robust resistance to overfitting . eoyeO ) = OcyeayeO=ya yuyea ycycyea . eoyeO ) . Furthermore, a Deep Neural Network (DNN) structure consists of an input layer, hidden layers, and an output layer, each of which comprises a set of neurons. Unlike Artificial Neural Networks (ANN), which typically have one or two hidden layers. DNNs feature a more significant number of hidden layers. This increased depth then enhances the network ability to generalize complex non-linear relationships between inputs and outputs . The output of a DNN is expressed in Eq. , where ycuycu represents the input, ycycn denotes the weights, ycaycn is the bias, and yce(O. is the activation function applied to the A key advantage of the DNN architecture lies in its capacity to capture and model intricate relationships between inputs and outputs . yeoyeO = OcycA yea=ya yeN. eoyeO y yeoyea yeEyeO ), yea OO . a, ycA] . This research employed the tree-parzen structured estimator (TPE) technique for hyperparameter tuning. Optimal hyperparameter determination is able to provide a better prediction model compared to default settings . The TPE technique is an enhancement of conventional Bayesian In conventional Bayesian methods, the surrogate function for determining hyperparameters uses a Gaussian function, whereas TPE employs probability density function (PDF) modelling to separate good and bad hyperparameter sets using kernel density estimation (KDE) . The process begins with the random initialization of hyperparameter combinations, which are then evaluated by an objective functionAiin this study, the root means square error (RMSE). Subsequently, the set of good and bad hyperparameters is determined using threshold estimation with KDE. A new configuration of the set of good hyperparameters is then found, and re-evaluated using the objective function. This process is performed iteratively n times . This research involved three to four hyperparameters In the SVR algorithm, three hyperparameters were optimized: regularization parameter (C), kernel, and kernel coefficient . For the kernel, a single option was used: the radial basis function (RBF). Meanwhile, for C and gamma, a normal distribution was used with a range of 0. 1 to 1000 for C, and 0. 0001 to 1 for gamma. In the DNN algorithm, there were four hyperparameters optimized: number of layers, number of neurons each layer, activation function, and The search space for the activation function was SoftMax and rectified linear unit . , as seen in Eq. and Eq. , where x is the input vector. yce ycuycn 0, ycu < 0 ycu, ycu Ou 0 cuycn ) = max (Ocycu yc=1 yce ycyceycoyc . cu ) = ycoycaycu ( ycuyc The search spaces for the DNN hyperparameter optimizer are adaptive moment estimation (Ada. and Adamax . variation of Adam for large parameter. There were four hyperparameters optimized in the GBR algorithm: the maximum number of trees for iteration . ax_estimator. , the maximum depth of each decision tree . ax_dept. , learning rate, and subsample. In the RF algorithm, four hyperparameters were optimized: the number of decision trees . ax_estimator. , the maximum depth of each decision tree . ax_dept. , the minimum number of samples for a data set to be split again . in_samples_spli. , and the maximum number of parameters used to build each DT . ax_feature. The max_features hyperparameter had three options: AusqrtAy, meaning that the number is the square root of the number of Aulog2Ay, meaning that the maximum number of features is log2. and None, meaning that all features are Table 2 presents the search spaces for each hyperparameter and algorithm. The prediction model was evaluated by comparing the predicted values with the actual ones and quantifying the results using five evaluation metrics, including correlation coefficient (CC), coefficient of determination (RA), root mean square error (RMSE), mean square error (MSE), and mean absolute error (MAE). The equations, value ranges, and optimal values for each metric are presented in Table 2, where ycycaycn represents the i-th actual data point, ycycyycn is the i-th predicted data point, yc iycy is the mean of the predicted data, and yc iyca is the mean of the actual data . Table 2. Model evaluation metrics Metric Formula Range Best value -1,1 Oc. cycyycn Oe ycycaycn )2 ycA OoycAycIya ycycyycn Oe ycycaycn | ycA Oc. cycaycn Oe yc ycyycn Oe yc ycy 2 Oc. c Oe yc cycaycn Oe yc ycy MSE OcycA ycn=1. cycyycn Oe ycycaycn ) ycA Ocycn=1. cycaycn Oe yc ycn=1 RMSE MAE ycA ycn=1 Cahyono et al. / Communications in Science and Technology 10. 135Ae147 Future epoch prediction The evaluated prediction model was then utilized to predict future epochs. Since the dataset generation process used data from the previous five epochs to predict a specific epoch, the future epoch prediction process also followed a stepwise approach using five epochs beyond the last observed data. The future prediction operated iteratively: the previous five epochs of data were used to predict one epoch ahead . cycycoycaycyc 1 ). Subsequently, to predict ycycycoycaycyc 2 , the predicted ycycycoycaycyc 1 was incorporated as ycuycOe1 . This process continued iteratively until ycycycoycaycyc 5 was obtained. intensity in the region. These results are consistent with previous research, which linked local rainfall phenomena with the global phenomena of El Niyo and La Niya . Results and Discussion Predicted data characteristics Precipitation prediction in this study involved thirteen parameters obtained from time series observation data. The prediction was naturally based on previous time series rainfall observations, allowing for the identification of its characteristics and behavioral patterns. The characteristics and behavior of rainfall can be analyzed in great detail based on historical daily precipitation data. In contrast, the analysis of global rainfall characteristics can be identified based on monthly and annual rainfall data. Fig. 4 shows the monthly precipitation characteristics from 2018 to 2023. The visualization shows that precipitation in the study area tended to be highest from December to March with the peak average occurred in February at 13. 75 A 8. 68 mm/day. Conversely, the lowest average precipitation was observed in July at 0. 79 A 2. mm/day. Fig. Characteristics of annual rainfall . ata source: CHIRPS accessed and processed using GEE, by the author. Fig. Monthly Precipitation Characteristics . ata source: CHIRPS accessed and processed using GEE, by the author. Fig. 5 shows that the overall distribution of rainfall values remained relatively consistent across the years, with only slight variations in rainfall trends between them. The annual rainfall distribution revealed that 2022 experienced the highest average rainfall at 7. 80 A 7. 92 mm/day and the highest maximum yearly rainfall at 48. 09 mm/day. Conversely, 2023 recorded the lowest average annual rainfall at 4. 79 A 6. 69 mm/day and the lowest maximum yearly rainfall at 31. 45 mm/day. These variations were associated to the ENSO phenomenon. In 2022, the ENSO index was predominantly positive, indicating a La Niya event that increased rainfall intensity in the equatorial region, including the study area. In contrast, the ENSO index in 2023 shifted negative, signifying an El Niyo event, reducing rainfall Fig. Prediction data variables The time-series behavior of each variable is depicted in Fig. To determine the relationship between each variable and rainfall in the study area, correlation analysis was Cahyono et al. / Communications in Science and Technology 10. 135Ae147 The results are presented as a correlation heatmap in Fig. The first row of the heatmap illustrates the correlation between rainfall values and other variables . Based on this visualization, all variables exhibited low to moderate correlations with rainfall with both positive and negative relationships observed. Of 12 variables, 8 showed a moderate correlation with rainfall including SST. WSU. CTH, and CBH with positive correlations, and WSV. SSR. CTP, and CBP with negative ones. In contrast, both LST . ay and nigh. and STMin/STMax had low correlations with rainfall in the study area with absolute correlation values below 0. process using TPE. Hence, 50 experiments were performed to find an optimal value for each hyperparameter. The TPE technique is able to reduce computational load as it only performs n-experiments to find the optimal hyperparameters, rather than trying each combination of hyperparameters one by Table 3. Hyperparameter tuning result using TPE each algorithm Hyperparameter Search spaces Optimum value SVR algorithm Regularization parameter (C) uniform . 1, 1. Kernel AorbfAo AorbfAo uniform . e4, . Kernel . DNN algorithm n_layers_dnn 6,7,8,9,10 n_unit_neurons 64, 128, 256, 512 activation_function AoreluAo. AosoftmaxAo AoreluAo AoAdamAo. AoAdamaxAo AoAdamAo Optimizer Fig. Correlation heatmap matrix among variables The characteristics of each data serve as the essential variable in the multivariate rainfall prediction modeling. Ideally, selecting variables with moderate correlations . oth positive and negativ. to rainfall can enhance the accuracy of the predictions . Based on the correlation analysis. WSU was found as the variable with the strongest positive correlation to rainfall and a correlation coefficient of 0. Conversely. WSV became the variable with the strongest negative correlation to rainfall with a correlation coefficient of -0. The results of correlation test indicated no significant relationship between rainfall pattern similarity and the spatial resolution of the data used. Of all components, temperatureAi including Land Surface Temperature (LSTD & LSTN) and Soil Temperature (STMin & STMa. Aishowed the weakest correlation with rainfall and most other components. contrast, cloud properties (CTH. CTP. CBH, and CBP) demonstrated a moderate correlation with wind speed . and v component. , shortwave radiation (SSR), sea surface temperature (SST), and precipitation. Results of multivariate rainfall prediction modeling Multivariate rainfall prediction was conducted using variables with an absolute correlation value greater than 0. Consequently, only nine variables were included: SST. WSU. WSV. SSR. CTH. CTP. CBH, and CBP. For each variable, data from the previous five-time steps . -1, t-2, t-3, t-4, and t-. were used as inputs for prediction. This resulted in a total of 45 predictive variables. The algorithms employed in this study were SVR. DNN. GBR, and RF, all of which were implemented with best hyperparameter architecture using TPE results. Table 3 shows the search spaces and optimal hyperparameters for each algorithm. The optimal hyperparameters were selected based on the combination that produced the minimum objective function value. In this study, there were 50 iterations for performing the optimization GBR algorithm n_estimators range . ,1000,. max_depth range . ,10,. Learning_rate uniform . 01, 0. 5, . RF algorithm n_estimators range . ,1000,. max_depth range . ,10,. min_samples_split range . ,10,. max_features AosqrtAo. Aolog2Ao. None None The models generated by each algorithm were applied to both training and testing datasets. Fig. 8 illustrates the application of each algorithm to the training data with a focus on rainfall data from January to March 2022 for more precise The results showed that most algorithms effectively captured the general pattern of the actual data during Of four algorithms. GBR provided the best visual fit, showing more minor discrepancies between predicted and actual values compared to others. The prediction results obtained through the DNN algorithm on the training data exhibited a relatively larger gap compared to actual data relative to other models. Other algorithms could consistently learn daily rainfall patterns effectively, particularly the GBR model. Most time steps were able to predict with great precision against the actual data. This condition was not entirely favorable as the prediction results on the testing data must be re-investigated. Since the gap with the actual data was very small, the GBR model could generalize the data well. The time-step segments in Fig. 9 showed that the DNN model tended to produce lower prediction values compared to the actual data with the most significant gap occurred during the first 15 days of March 2023. Fig. 9 shows the results of applying each algorithm's prediction model to the testing data. Similar to the training data visualization, the comparison between actual and predicted values was focused on the time period from January to March The visualization indicated that all prediction models Cahyono et al. / Communications in Science and Technology 10. 135Ae147 produced patterns closely aligned with the actual data, suggesting that they are generally reliable when applied to new In the time-step segment, the SVR model was found consistent in predicting actual data with greater precision than other models. While the prediction values GBR model were nearly aligned with actual values in the training data visualization, it showed a significant gap compared to the SVR model in the testing data. Visually, the GBR model was found less consistent in accurately predicting actual values. The DNN model exhibited a similar pattern in the testing data with a relatively larger gap compared to other prediction models. Prediction values in the testing data tended to be lower than actual values. The difference between actual and predicted values was quantified by the residuals, representing a difference between the actual and predicted data. Fig. 10 illustrates the residual distribution for each algorithm on the training data. Of the algorithms, the SVR algorithm had the smallest median value 116y10-3, while RF had the highest one at -0. Although SVR had the smallest median value. GBR had the smallest standard deviation with a value of A0. 41, meaning it had the most stable relative bias. This can be seen in how the outlier points in the GBR model clustered near the interquartile range (IQR) and formed a shorter boxplot compared to other The median value in the GBR model was -0. DNN was the model with the highest median and standard deviation values with the values of 1. 726 and A2. 492, respectively. This can be visually seen where the DNN model formed a relatively long boxplot compared to other models. The residual values of each algorithm for the testing data are illustrated using a box plot in Fig. Of the models. SVR remained the model with the lowest median residual value of In addition, it was the one with the lowest residual standard deviation of A0. The boxplot generated by the SVR model was visually shorter than those generated by other In contrast, the DNN model had the highest median and standard deviation values with the values of 0. 950 and A2. 132, respectively. This model visually formed the longest boxplot of other models. Although the GBR model had better residual stability on the training data, when compared to the testing data. SVR outperformed the residual stability. Fig. Prediction results of each algorithm on the training data Fig. Residuals of each algorithm on the training data Fig. Residuals of each algorithm on the testing data Prediction model evaluation Fig. Prediction results of each algorithm on the testing data The prediction model was evaluated using both training and testing data, employing five evaluation metrics: CC. RA. RMSE. Cahyono et al. / Communications in Science and Technology 10. 135Ae147 MAE, and MSE. The results for each algorithm and dataset are summarized in Error! Reference source not found. On the t raining data, the SVR algorithm achieved the highest CC and RA values, 0. 982 and 0. 965, respectively, indicating its strong correlation and goodness-of-fit. Conversely, the DNN algorithm recorded the lowest ones, 0. 943 and 0. For error-based metrics (RMSE. MAE, and MSE). GBR demonstrated superior performance with the lowest values at 0. 426, 0. 332, and 0. 182, respectively. DNN exhibited the highest ones at 3. 363, 2. 351, and 312, respectively. well than the ML algorithm. It did not adequately learn the daily rainfall fluctuations by considering multiple variables. However, when using a very large dataset the DL algorithm can perform well . Therefore, a longer training data series is needed to learn the rainfall fluctuations in the study area. Table 4. Evaluation Metrics for Each Algorithm Metric RMSE MAE MSE Dataset SVR DNN GBR Training Test Training Test Training Test Training Test Training Test The SVR algorithm performed better when evaluated with testing data. It had the highest CC and RA values of 0. 974 and 948, respectively but showed the lowest error-based metrics with RMSE. MAE, and MSE values of 1. 366, 0. 947, and 1. While, the DNN algorithm consistently demonstrated the lowest performance with CC and RA values of 932 and 0. 890, respectively. The RMSE. MAE, and MSE values for the DNN model were 2. 804, 2. 351, and 7. These results highlighted the variability in algorithm performance between training and testing datasets . The findings of this research indicated that the SVR algorithm was the best algorithm compared to the DNN. GBR, and RF algorithms. It consistently produced prediction values very close to the actual data. A study by Nayak et al. produced the RMSE values of 0. 912 to 1. 091 for monthly rainfall predictions using the SVR algorithm . In contrast. Wang et al. showed that SVR performed worse than RF and Bayesian ridge regression (BRR) for univariate rainfall prediction . The GBR algorithm outperformed the SVR algorithm on the training data, but did not consistently outperform it on the testing data evaluation. Similar to the study by Sumith . GBR had an almost perfect RA of 0. 95 on the training data, but when evaluated on the testing data it had an RA value of 0. Nevertheless, in this study the difference in evaluation results between the training and testing data for the GBR model was less significant with an RA value on the testing data greater than 0. It was because the GBR algorithm continuously iterated residual values, allowing the model to closely approximate actual values in the testing data. consequence, the model performed poorly in generalizing unseen data . The DL algorithm in this study performed less Fig. Box-plot of KFCV results evaluated by MAE Fig. Box-plot of KFCV results evaluated by R2 This research employed k-fold cross validation (KFCV) to evaluate the stability and robustness of each prediction model. KFCV was conducted by setting the number of folds to 10, thereby dividing the entire dataset into 10 folds of equal size, with 1-fold designated as the testing data and the remaining 9 folds as the training data. The fold used as testing data was used sequentially, and the model was consistently fitted to the training fold in combination with other nine folds. The results of the KFCV process are visualized in the form of box plots, as shown in Fig. 12 and Fig. The metrics used for KFCV were MAE and RA. In the box-plot visualizations, the x-axis represents the prediction model, and the y-axis represents the metric value, with each model grouped into training and testing The training dataset was visualized with a blue box, and the testing data was shown in red. The whiskers represented the range of minimum and maximum metric values. The results of KFCV showed that the MAE value of the GBR model and training data tended to have the smallest MAE value compared to other models and exhibited high stability, as indicated by the short whiskers produced. The MAE value of the GBR model's Cahyono et al. / Communications in Science and Technology 10. 135Ae147 training data was 0. 35A0. Meanwhile, when all models were evaluated using the testing data, the SVR model was found more robust, producing the smallest average value and relatively higher stability compared to other models. The short whiskers in the boxplot represented these characteristics. The average MAE value for the SVR model was 1. 14 A 0. The KFCV R2 metric results followed the same pattern as the MAE. The GBR model had the highest value with high robustness at 00 A 0. Meanwhile, the SVR model had the highest R2 and robust on the testing data with a value of 0. 94 A 0. Therefore, in the present study, the SVR model demonstrated both the best performance and strong robustness across various splitting The performance of each prediction model was further evaluated using a Taylor diagram, visually representing the relationship and comparison between prediction and actual A Taylor diagram includes three main components: the radial axis, which indicates the standard deviation of the data. isolines, which denote the RMSE values. and the quarter-circle arcs, which represent the CC. The reference point on the diagram corresponds to the actual values, serving as the benchmark for assessing the closeness of the predicted values to the actual data. The Taylor diagrams for the training and testing data are illustrated in Fig. 14 and Fig. 15, respectively. Fig. 14 shows that all algorithms exhibited a strong correlation with the actual data on the training dataset (CC>0. The GBR algorithm had the closest distance to the actual standard deviation value, the lowest RMSE value, and the highest correlation. In contrast, the DNN algorithm had the opposite characteristics, namely the furthest distance from the actual standard deviation value, the highest RMSE, and the lowest correlation. The Taylor diagram for the testing data (Fig. reveals that SVR not only maintained low RMSE and high CC values but also aligned closely with the standard deviation reference line of the actual data. Based on both quantitative and visual evaluations, the SVR algorithm outperformed other algorithms in this study, showcasing its robust and reliable predictive This research explored the potential for overfitting using learning curves. A learning curve is a graph that illustrates the relationship between the percentage of total training data used and the resulting accuracy, both on new training data . aken from 0-100% of the training dat. and testing data. The testing data used to evaluate each scenario was the testing data resulting from an 80:20 split in the initial process. Meanwhile, the fraction of the training data was determined based on the 80:20 split during the initial process. In each scenario, the model was fitted with the resulting training data fraction. Fig. 16 presents the learning curve analysis results for each model that was visualized by different colors. Meanwhile, the difference between the training and testing data results for each algorithm was visualized by the different types of lines, where training data and testing data were represented by a solid line and a dashed line, respectively. The metric used for learning curve analysis was R2. The results obtained then showed that the more training data used . he data fraction approaching . , the more convergent the graph between training and testing data for each model. This indicated no significant gap between the training and testing data. However, when compared relatively across models, the RF and GBR algorithms showed the highest gap between training and testing data compared to SVR and DNN, indicating the potential of these two models for Nevertheless, the R2 value for the testing data of both algorithms was above 0. 8, meaning the models are still able to make good predictions on unseen data. Fig. Taylor diagram for training data Fig. Taylor diagram for testing data Fig. Learning curve each model Cahyono et al. / Communications in Science and Technology 10. 135Ae147 Future epoch prediction Subsequently, the prediction model was applied to forecast rainfall values for future epochs. Five future epochs were set to align with the number of lags used in the dataset generation This corresponds to predicting rainfall for the next five Future epoch forecasts were made recursively, where the output at epoch t 1 became the input for the forecast at epoch t 2, and thus onward. Fig. 17 illustrates the forecasting results for these five future epochs, as produced by each algorithm. The data were visualized in a line chart accompanied by a box plot to illustrate the distribution of data predictions across obtaining an RA of 0. Overall, this study highlighted SVR as the best-performing algorithm and the one that was resistant to overfitting. The limitation of this research is that it merely provided predictions in the time domain, rather than spatial domain predictions. Nonetheless, the results are promising for daily rainfall predictions in the study area, utilizing open-access remote sensing big data. Despite differences in the spatial resolution of the data, the accuracy of the predictions is It is essential that the data have the similar temporal resolution . This research is expected to apply across various sectors, particularly hydro-meteorological disaster management. Acknowledgements Fig. Future forecasting for each epoch Predictions were generated for the period from December 31, 2023, to January 5, 2024. During epoch t 1, the four algorithms predicted rainfall closely, as indicated by the small box plots generated. The standard deviation value at epoch t 1 was A1. However, during epochs t 2 to t 5, the DNN algorithm had prediction values relatively far from the SVR. RF, and GBR prediction values. A significant difference was observed at epoch t 2, with a standard deviation of A4. The t 2 epoch had the highest average predicted rainfall at 14. While, the lowest rainfall was recorded at t 5 with an average of 3. 467 mm. We would like to extend our heartfelt acknowledgments to the research team of "Smart Sustainable Urban Flood Management: A Decision Support System Towards a Sensing City Ae Lessons Learned from the Surabaya Case for IKN Development", for their collaboration on this joint research. Our deepest appreciation also goes to the Faculty of Engineering. Universitas Gadjah Mada (UGM), for providing research funding, facilities, and various forms of support that have facilitated the successful progress of this research. References