Modernizing harvested paddy areas estimation through satellite imagery time series and machine learning

2.1. Satellite imagery time series

This analysis utilizes Satellite Imagery Time Series (SITS) data derived from historical Sentinel-1 imagery, obtained via the Alaska Satellite Facility (ASF) API for the period 2021–2023, at a biweekly temporal resolution. ⁹ Instead of the satellite image in one time, the longitudinal data of satellite image, which is SITS, is chosen to enhance the modelling of temporal trends. ¹⁰ To ensure each image quality, standard preprocessing steps—such as orbit file correction, thermal noise and border noise removal, radiometric calibration, speckle filtering, terrain correction, and mosaicking—is applied.^11,12 All preprocessing procedures are conducted on an in-house server using the ESA-SNAP Python API, producing a time series of Sentinel-1 imagery focused on the VH-band, with a 20-meter spatial resolution. ¹¹ After preprocessing, integration between Sentinel-1 imagery data and ASF field observations for the period February 2022–December 2023 is performed by spatial join techniques. Field observations are collected from 300 m × 300 m plots, each subdivided into nine 100 m × 100 m sub-grids. The corresponding time series segments are illustrated in Figure 1, where the yellow shaded region represents the final 12 weeks (± 90 days), corresponding to the most recent paddy cultivation period.

OPEN IN VIEWER

Figure 1.

The ASF plot and sentinel 1 integration.

The ASF data, collected from monthly observations, is used as the ground truth label for the classification task. Phenological stages are recorded at the centroid of each sub-grid using geotagged points. Each monthly observation is matched (paired) with a set of SITS data of length 31, corresponding to the total images of Sentinel 1 collected biweekly in 1 years. However, several adjustments are made to stage labeling: non-paddy and no-vegetation categories are combined, and consecutive observations of land preparation and harvest are recategorized, as detailed in Table 1. For consecutive land preparation, only the last observation is designated as “land preparation,” while others are assigned to a newly created class (“bare-land”). Conversely, for consecutive harvest records, only the first observation is categorized as “harvest,” with the remainder also classified as “bare-land.”

3.1. Preliminary

As a preliminary step, represented in Figure 2, the regencies in focus regions undergo a clustering process based on agricultural and terrain characteristics, involving data as follows: previous total harvested area, previous total production area, number of paddy plantation cycles, elevation, evapotranspiration, and Land Surface Temperature (LST). The clustering process is designed to generate distinct groups of regencies that serve as the basis for machine learning modelling respective to the spatial constraints. Consequently, rather than employing a single global model, multiple models are developed, each tailored to the characteristics of its respective cluster. For this purpose, the Spatial ‘K'luster Analysis by Tree Edge Removal (SKATER) algorithm is implemented due to its competence in handling spatial clustering. ¹³

OPEN IN VIEWER

Figure 2.

Preliminary analysis step.

After the first step is applied, three main preprocessing procedures are undertaken: imputing missing data, heterogeneity filtering, and noise filtering to ensure a well-representative training set.

The preprocessing step, depicted in Figure 3, begins by addressing missing values in the Sentinel data using the Whittaker-Eilers smoothing algorithm, which effectively fills gaps in time series data. ¹⁴ This ensures data continuity and reliability before further processing. Subsequently, to improve training dataset quality, grids with high heterogeneity are excluded using variance-based filtering. Sub-grids with excessive intra-cluster variability, measured by the variance of VH values over time, are removed if they exceed the third quartile (Q3) within each cluster. This criterion reflects mixed cropping practices commonly observed in Indonesia, which result in elevated VH variability within individual sub-grids. While this filtering step removes a non-negligible proportion of observations, it represents a deliberate trade-off between sample size and data homogeneity. By prioritizing more homogeneous agricultural patterns, the resulting training dataset improves model stability and robustness, at the cost of reduced representation of highly heterogeneous areas.

OPEN IN VIEWER

Figure 3.

Data preprocessing step.

Lastly, Self-Organizing Map (SOM) filtering ¹⁵ is applied to noise filtering and enhance data representativeness. Each cluster is assigned 100 neurons, and neurons are selected based on the majority phenological stage. Observations most similar to these selected neurons are retained, while others are discarded.^16,17 To introduce controlled randomness and reduce overfitting, 5% of the discarded observations are reintroduced into the training set.

Model development is conducted to several process, as depicted in Figure 4. The modelling process is applied in each regional cluster to improve the accuracy. At the first stage, the dataset is divided into two subsets: one containing all observations and another with 70% of the data. The smaller subset is then used to train the machine learning model, which is subsequently applied to predict the full dataset. This subset is further split into training and validation sets with a 70:30, while ensuring the distribution of paddy phenological stages remains representative. After thorough preprocessing of the training dataset, the modelling phase begins with hyperparameter tuning and cross-validation to optimize performance. ¹⁷ XGBoost is selected as the algorithm due to its robustness, efficiency, and ability to capture complex patterns. Since the task focuses on multiclass classification with imbalanced categories caused by the uneven distribution of paddy phenological stages, the imbalance is preserved to reflect real-world conditions.

OPEN IN VIEWER

Figure 4.

Model development step.

Hyperparameter tuning is conducted using the Optuna library, which efficiently explores parameter ranges to optimize model accuracy. ¹⁸ The macro-averaged F1 score is chosen as the optimization objective since it is well-suited for imbalanced classification tasks. ¹⁹ A three-fold cross-validation is applied to ensure stability, with tuning centred around critical hyperparameters such as ‘max_depth’, ‘learning_rate’, ‘reg_lambda’, and ‘min_child_weight’. The model uses ‘multi:softprob’ as the objective function for multiclass probability outputs, and early stopping (50 iterations) is implemented to reduce overfitting, supplemented by weighted sampling strategies for handling imbalance (Chen & Guestrin, 2016). Once the optimal parameter set is identified, it is used to finalize the model through a 5-fold cross-validation setup for robust evaluation. To further prevent overfitting, boosting iterations are capped at 10,000 with an early stopping threshold of 1000 rounds. Finally, performance metrics including overall accuracy, F1-Macro, F1-Micro, ROC-AUC, and phase-specific accuracies are averaged across folds to assess generalization.

Using the models developed for each cluster, predictions are generated for all observations from February 2022 to December 2023, with the process depicted in Figure 5. These predictions focus solely on the centroids of each sub-grid, ensuring consistency with ASF observation data, which is also centred on centroids. Additionally, predictions are made at locations corresponding to ASF sample points, effectively simulating the ASF data collection process through machine learning predictions based on ground-truth data. This alignment minimizes discrepancies between predicted and observed data.

OPEN IN VIEWER

Figure 5.

Prediction and estimation step.

The predicted data is then used to estimate paddy phenological stages, following ASF's dot sampling procedure. This approach ensures consistency with ASF's estimation methods. Both point estimates and sampling errors are calculated to evaluate prediction accuracy. The dot sampling method maintains alignment with ASF's approach, while calculating sampling errors provides insights into prediction reliability. Confidence intervals from ASF and this study are compared to assess the level of agreement between the two methods. High alignment suggests that the machine learning model effectively replicates observed data and maintains consistency with field-based measurements, demonstrating the model's reliability and accuracy.

4.1. Clustering results

Figure 6 presents the clustering results derived from its agricultural characteristics. Cluster 0, located northwest, is characterized by low elevation, high land surface temperature (LST), low harvested area, and low production, indicating urban and coastal environments. Cluster 1, located on the north, features low elevation and moderate land surface temperature (LST), with intensive paddy cultivation and high productivity, signifying a strong agricultural presence. Cluster 2, located in the northeast, exhibits both urban and agricultural characteristics, with high land surface temperature (LST), and moderate paddy production. Cluster 3, located in the southeast, has moderate elevation but high evapotranspiration, indicative of dense vegetation and less productive agriculture. Cluster 4, located in the south, is characterized by high elevation and low land surface temperature (LST), suggesting fertile highlands yet paddy plantation is not the primary sector. Lastly, Cluster 5, located in the centre, shares similar highland characteristics with abundant vegetation but limited paddy farming.

OPEN IN VIEWER

Figure 6.

Clustering results in West Java, Indonesia.

Following the cluster results, the data collected undergo the modelling after preprocessing is conducted and Figure 7 represents the data distribution on the training dataset. As depicted in the figure, the presence of imbalanced data can be noted. Cluster 1 consistently exhibits the highest counts across most observation classes excluding the non-paddy and non-vegetation classes. In contrast, Cluster 3 consistently has the lowest counts, primarily for Vegetative 2 and the land preparation classes. The other cluster displays moderate variations for the observation classes, although the highest count of non-paddy and no-vegetation classes is shown in both clusters 0 and 5. The data itself is further modelled with results as presented in Table 2.

OPEN IN VIEWER

Figure 7.

Data distribution on training data.

OPEN IN VIEWER

Table 2.

Accuracy metrics on training and validation dataset.

Training Dataset	Cluster	Accuracy	F1 Macro	F1 Micro	ROC AUC	Accuracy [Generative]	Accuracy [Paddy Phase]
	0	0.9201	0.9015	0.9201	0.9855	0.9454	0.9519
	1	0.8958	0.8794	0.8958	0.9842	0.9136	0.9211
	2	0.9319	0.9221	0.9319	0.9912	0.9523	0.9532
	3	0.6028	0.5729	0.6028	0.8761	0.9747	0.7451
	4	0.8806	0.8425	0.8806	0.9687	0.957	0.896
	5	0.8096	0.7456	0.8097	0.9474	0.9644	0.8777
Validation Dataset	Cluster	Accuracy	F1 Macro	F1 Micro	ROC AUC	Accuracy [Generative]	Accuracy [Paddy Phase]
	0	0.9281	0.9113	0.9281	0.9879	0.9494	0.958
	1	0.8966	0.8881	0.8967	0.9857	0.9091	0.9183
	2	0.9315	0.922	0.9315	0.9911	0.9486	0.9526
	3	0.6109	0.5849	0.6109	0.8842	0.9746	0.7455
	4	0.8867	0.8515	0.8867	0.9726	0.9616	0.9031
	5	0.8096	0.7456	0.8097	0.9474	0.9644	0.8777

In the training data results, Clusters 0, 1, 2, and 4 achieved accuracy and F1 macro scores exceeding 85%, indicating that the XGBoost model performs optimally in these regions regardless of class imbalance between paddy growth phases. Furthermore, the AUC values for these clusters surpassed 98%, demonstrating the model's strong ability to distinguish between classes. The accuracy for both the generative phase and paddy phase also remained above 85% across these clusters, reflecting the model's consistent performance in identifying paddy growth phases. Conversely, Cluster 3 exhibited significantly lower performance, with an accuracy of 60.28% and an F1 macro score of 57.25%. This suggests that additional attention is required to enhance the model's predictive capability in this region. Cluster 5 also showed suboptimal results, with accuracy and F1 macro scores below 85%, despite achieving a generative phase accuracy of 96.44%. This indicates that while the model struggles overall in this cluster, there is potential for improvement, particularly in the generative phase.

In the validation data results, Clusters 0, 1, 2, and 4 again demonstrated the highest performance, with accuracy and F1 macro scores exceeding 85%. Their AUC values remained above 98%, reaffirming the model's strong classification capability. The accuracy for both the generative phase and the paddy phase was also consistently high, exceeding 90%. However, Cluster 3 continued to show the lowest performance, with an accuracy of 61.09% and an F1 macro score of 58.49%, indicating persistent difficulties in recognizing patterns in this region. Similarly, Cluster 5 maintained accuracy and F1 macro scores below 85%, aligning with the training data results, thus requiring further refinement to enhance model performance, particularly in the paddy phase.

Based on these findings, Clusters 0, 1, 2, and 4 emerge as the most suitable regions for model implementation due to their consistently high performance across all metrics. In contrast, Cluster 5, while underperforming, exhibits potential for improvement, especially in the generative phase, suggesting that further tuning could enhance its predictive capability. These results highlight the necessity of targeted model adjustments to optimize classification accuracy across all clusters, particularly in regions with lower performance.

The model's performance aligns closely with the geographical characteristics of each cluster. Clusters 0, 1, 2, and 4, which are primarily flat or low-elevation regions, demonstrated high accuracy and F1 macro scores exceeding 85%, indicating the model's effectiveness in areas with distinct agricultural, mainly for paddy crops. Their AUC values above 98% further reinforce the model's reliability in distinguishing the paddy phenological stage. In contrast, Clusters 3 and 5, which are characterized by higher elevation, dense vegetation, and complex terrain, exhibited lower performance, with Cluster 3 performing the worst. This suggests that the model struggles in mountainous regions, likely due to environmental variability affecting paddy growth patterns.

Figure 8 depicts the distribution of the paddy phenological stage from February 2022 to December 2023, measured in hectares. Throughout the timeline, fluctuations in each phase are observed. Bareland and Land Preparation exhibit seasonal peaks, notably increasing towards the end of each year. The Vegetative 1 and Vegetative 2 phases follow a similar trend, indicating the progression of paddy growth, with peaks preceding the Generative phase. The Generative phase shows moderate variation, aligning with prior vegetative growth. Harvested areas fluctuate in response to these phases, peaking after the Generative phase. The Non-Paddy and Non-Vegetation categories remain relatively stable.

OPEN IN VIEWER

Figure 8.

Paddy's phenological stage derived from the prediction.

Figure 9 presents the estimated harvested area derived from the machine learning model and its comparison with ASF data from February 2022 to December 2023. The model reproduces the overall temporal and seasonal patterns observed in the ASF estimates, indicating a strong ability to capture the main dynamics of harvested area. In particular, the primary harvesting peaks in March 2022 and March 2023 are clearly identified, as well as the decline after October 2022 and the inflection point in December 2022 preceding the subsequent increase toward the 2023 peak. The predicted series also reflects secondary peaks in July of both years, which broadly align with fluctuations observed in the ASF data.

OPEN IN VIEWER

Figure 9.

Estimation with confidence interval of harvested area from Feb 2022–Dec 2023.

However, the model exhibits systematic overestimation during several periods, most notably in September 2022, where a pronounced increase is predicted that is not supported by the ASF estimates. As a result, the confidence intervals of the model-based estimates do not overlap with those of the ASF data for substantial portions of the study period. Although overlapping confidence intervals are observed in some months, suggesting partial agreement between the machine-learning estimates and ASF data, this agreement is not consistent over time. Extended periods of overestimation lead to non-overlapping confidence intervals, indicating that the discrepancies between the two approaches are better characterized as moderate rather than minor.

Therefore, the agreement between the machine learning estimates and ASF data should be interpreted primarily in terms of their shared seasonal behavior, rather than as full statistical equivalence. Differences in the width of the confidence intervals further indicate varying levels of uncertainty, with ASF estimates exhibiting greater variability in several months.

Turning to Figure 10, the relative standard error (RSE) analysis provides insight into the precision of both ASF and machine learning–based estimates. The RSE values for both sources remain below 10%, indicating that the estimates meet commonly used quality thresholds for statistical reporting. This threshold is consistent with United Nations guidelines, which recommend relative error margins between 5% and 10% for key indicators intended for analytical and policy applications. ²⁰

OPEN IN VIEWER

Figure 10.

Relative standard error (RSE) of harvested area from feb 2022-Dec 2023.

The machine learning–based predictions generally exhibit lower RSE values than the ASF estimates, indicating a higher degree of internal consistency and temporal stability. However, lower RSE values should not be interpreted as higher accuracy. In several periods, the model produces systematically overestimated values with relatively high certainty, reflecting reduced variability rather than reduced bias. Consequently, while the model delivers stable and precise estimates, these estimates may still deviate from field-based measurements.

Overall, the results highlight the model's strengths in generating stable, low-cost, and temporally consistent estimates that effectively capture seasonal harvesting patterns. At the same time, the presence of systematic bias underscores the need for cautious interpretation when using such outputs for official statistics. Machine learning–based estimates are therefore best viewed as a complementary information source, rather than a direct replacement for field-based data, unless appropriate bias correction or calibration procedures are applied.