# Reconstructing satellite temporal series data under cloudy conditions: Application in predicting rice growth phases

B Suseno, Guilhem Brunel, H Wijayanto, K Sadik, F M Afendi, Bruno Tisseyre

<figure><img src="image_1.png" /></figure>

## To cite this version:

B Suseno, Guilhem Brunel, H Wijayanto, K Sadik, F M Afendi, et al., Reconstructing satellite temporal series data under cloudy conditions: Application in predicting rice growth phases. Smart Agricultural Technology, 2025, 12, pp.101378. 10.1016/j.atech.2025.101378. hal-05280364

HAL Id: hal-05280364

https://hal.inrae.fr/hal-05280364v1

Submitted on 24 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Reconstructing satellite temporal series data under cloudy conditions: Application in predicting rice growth phases

B. Suseno<sup>a, b</sup>, G. Brunel<sup>b</sup>, H. Wijayanto<sup>a, *,</sup>, K. Sadik<sup>a</sup>, F.M. Afendi<sup>a</sup>, B. Tisseyre<sup>b, c</sup>

## ARTICLE INFO

### Keywords:

- Remote sensing

- Crop phenology

- Spatiotemporal reconstruction

- Multistep smoothing

- Interval-based approach

## ABSTRACT

Vegetation indices derived from satellite observations, particularly the NDVI, are widely used for monitoring crop development and estimating yields. In tropical regions such as Indonesia, however, persistent cloud cover causes large gaps, undermining the reliability of NDVI time series and limiting their agricultural utility.

This study proposes a two-stage spatiotemporal reconstruction framework to enhance NDVI time series under such conditions. The approach combines spatiotemporal filling with multistep temporal smoothing through an interval-based strategy that accommodates diverse planting schedules and performs effectively on shorter time series.

Experiments with simulated datasets demonstrated that the method consistently outperformed benchmark techniques, achieving the lowest reconstruction errors compared to the Savitzky-Golay filter and cubic spline smoothing. Applied to Sentinel-2 NDVI from 2,460 rice field sites across three Indonesian regions, the method enhanced rice growth phase classification with a random forest model, improving precision by up to 19%, recall by 24%, and F1-score by 21%. These results demonstrate the framework's robustness in handling large cloud-induced gaps and its potential for near-real-time agricultural monitoring.

## 1. Introduction

Rice plays a strategic role in Indonesia's food security, making accurate monitoring of its growth phases essential for yield forecasting and agricultural policy. The Area Sampling Frame (ASF) survey remains an important field-based monitoring approach [33], but it is labor-intensive and often impractical in inaccessible or rugged areas [13,23]. Satellite remote sensing offers a compelling alternative due to its capacity for consistent, large-scale observations of land surface dynamics [40]. In particular, time-series satellite imagery enables phenological analysis and improves monitoring efficiency in remote agricultural regions [10,22,36].

Vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), derived from multispectral satellite sensors, are widely used to monitor vegetation greenness [41]. Sentinel-2A and 2B, with their spatial resolutions of 10–60 m and a revisit cycle of five days, are especially well suited for monitoring agricultural dynamics [11,12], However, optical sensors are highly susceptible to cloud cover, particularly in tropical regions, resulting in large gaps in the data record [15, 42]. These gaps impair crop phenology assessments, especially in

regions with prolonged cloudiness [4,14,23,24]. Consequently, robust gap-filling or reconstruction methods are essential to ensure reliable satellite-based crop monitoring [28,39].

Numerous reconstruction techniques have been developed, broadly categorized into temporal-based, frequency-based, and hybrid approaches [25]. Temporal methods such as the Savitzky-Golay filter smooth each pixel independently but struggle with large gaps, as they rely on local temporal continuity [28,34,35]. This weakness is particularly problematic in tropical regions where cloud cover can persist for weeks [7,9,21,25]. Frequency-domain methods, such as Harmonic Analysis of Time Series (HANTS), model NDVI as a combination of Fourier components [31]. While effective in regular seasonal contexts, their performance declines in settings with asymmetric cycles or highly variable planting schedules, as commonly observed in Indonesian rice cultivation [1,3,8,19].

Hybrid approaches that combine spatial and temporal information have shown strong potential. Notable examples include Spatio-Temporal Tensor Completion [9] and the Spatio-Temporal Savitzky-Golay (STSG) method [7,25]. STSG, originally developed for MODIS NDVI, can effectively handle long missing sequences. Nevertheless,

<footnote><sup>a</sup> School of Data Science, Mathematics, and Informatics, IPB University, Meranti Wing Street, Dramaga, Bogor, West Java 16680, Indonesia</footnote>

<footnote><sup>b</sup> ITAP, Univ Montpellier, INRAE, Institute Agro, 2 Place Pierre Viala, Montpellier 34060, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: Bayu.suseno@apps.ipb.ac.id (B. Suseno), hari@apps.ipb.ac.id (H. Wijayanto).</footnote>
STSG faces key limitations. First, it requires multi-year NDVI data, as demonstrated in the original study, which used 16 years of data and required at least 15 days of processing on a single processor, rendering it impractical for large-scale operational applications. Second, it assumes that all pixels within a neighborhood follow similar seasonal NDVI trajectories. This assumption does not hold in agricultural systems like those in Indonesia, where staggered planting schedules create multiple cropping cycles within a single year. Third, it relies on the availability of cloud-free neighboring pixels, which are often unavailable in tropical regions during the wet season.

To overcome these limitations, a novel hybrid reconstruction method is proposed that integrates spatiotemporal pattern matching with multistep temporal smoothing. This study makes three key contributions: (i) a spatiotemporal filling strategy is introduced to reconstruct large cloud-induced gaps by leveraging historical temporal patterns and spatiotemporal information from other sites through an interval-based approach, ensuring robustness under persistent cloud cover and staggered planting schedules; (ii) a multistep smoothing procedure is developed to progressively refine the reconstructed series by aligning it with clearer observations, thereby preserving temporal continuity; and (iii) the practical utility of the framework is demonstrated through experiments on both synthetic datasets and real Sentinel-2 NDVI series from Indonesian rice fields. The rice growth phase was validated against benchmark reconstruction methods using field-based observations from the ASF survey.

## 2. Materials and methodology

### 2.1. Study area

This study was conducted across three Indonesian regions with contrasting agro-climatic conditions: Indramayu, Garut, and Pontianak (Fig. 1a), covering the period from January 2021 to December 2022. In total, 2,460 observation sites were selected. Indramayu (1,522 sites; Fig. 1c) represents Indonesia's largest rice-producing area, characterized by extensive irrigated lowland rice fields under a tropical monsoon climate. Garut (907 sites; Fig. 1d) is a highland region with predominantly flat terrain at higher altitudes. Pontianak (31 sites; Fig. 1b), while not a major rice-producing region, was included as a representative equatorial area, characterized by high solar radiation and biomass accumulation, providing a contrasting setting for evaluating the reconstruction method.

### 2.2. Data sources

#### 2.2.1. Rice growth phases

Monthly rice growth-phase observations were collected by trained surveyors as part of the ASF survey conducted by BPS-Statistics Indonesia. The ASF survey design is based on a regular 100 m × 100 m grid, with observation sites located at the centers of sampled grid cells. In total, 2,460 sites were visited, and their geographical coordinates were recorded during ground checks. These field data provided reference information on rice growth phases, which were classified into five

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Study area comprising three distinct regions in Indonesia: (a) Geographic locations of the study regions; (b) Observation sites in Pontianak, a lowland equatorial area characterized by high biomass; (c) Observation sites in Indramayu, Indonesia's largest irrigated rice-producing region; and (d) Observation sites in Garut, a highland area with predominantly flat terrain.</figcaption></figure>
categories (Fig. 2):

1. Vegetative 1 (early vegetative): Seed germination, emergence, and early tillering, when plants undergo rapid vertical and foliar growth but reproductive development has not yet begun.

2. Vegetative 2 (late vegetative): Continued tillering and stem elongation, during which the canopy typically reaches peak density and greenness.

3. Generative (reproductive): Panicle initiation, flowering (anthesis), and grain filling, marked by the emergence of the panicle and fertilization.

4. Harvest: Grain maturation and senescence, characterized by a shift in grain colour from green to golden/brown as maximum dry weight is reached.

5. Land Preparation Phase: Post-harvest operations such as ploughing, levelling, tilling, and puddling to prepare the soil for the next planting cycle.

### 2.2.2. Satellite images

NDVI and CloudScore+ values [29] were retrieved from Sentinel-2 imagery via the Google Earth Engine (GEE) platform using the coordinates of the 2,460 ASF sites. For each site, the pixel corresponding to the center of the 100 m × 100 m grid cell was extracted to ensure consistency with the ASF design. In total, 824,526 pixel-level observations were processed: Indramayu (427,441), Garut (351,985), and Pontianak (45,100). NDVI and CloudScore+ values were averaged for sites spanning multiple pixels to ensure consistency with field observations.

The NDVI was calculated from Sentinel-2 imagery using the red (RED) and near-infrared (NIR) bands of the Multispectral Instrument (MSI) sensor [32] (Eq. 1):

$$NDVI = \frac{NIR - RED}{NIR + RED} \qquad (1)$$

NDVI values range from -1.0 to 1.0 and serve as a proxy for vegetation greenness [41].

Pixel quality is critical when working with optical satellite data, which is frequently affected by cloud and shadow contamination. CloudScore+ is a detection metric that provides more accurate assessments than methods such as s2cloudless or CNN-based L3 [30]. It is based on a weakly supervised learning framework that uses an atmospheric similarity metric to bootstrap a globally distributed training dataset. This approach avoids reliance on fully annotated labels and enables the training of a space-time context model on millions of Sentinel-2 video clips [29]. CloudScore+ values range from 0 to 1, with

higher values indicating greater confidence in cloud- and shadow-free observations.

### 2.3. Influence of clouds and shadows on NDVI data

Indonesia's tropical climate is characterized by dry and wet seasons. The wet season, typically from September to February, is marked by persistent cloud cover, which frequently disrupts NDVI data by creating prolonged and recurrent data gaps. A CloudScore+ threshold of ≤ 0.1 was applied across all sites to identify severely affected observations.

Using this criterion, more than 40% of NDVI data were classified as severely contaminated during the study period, with proportions rising sharply in the wet season and reaching 88% in February 2022 (Fig. 3). Such disruptions pose major challenges for monitoring rice phenology, particularly during key growth transitions.

Fig. 4 further illustrates these effects. In Garut (Fig. 4a) and Indramayu (Fig. 4b), NDVI values on less cloudy days show clear seasonal dynamics, peaking during the Vegetative 2 phase (maximum canopy density) and declining to their lowest during the Land Preparation phase (minimal cover). In contrast, Pontianak (Fig. 4c) exhibits weaker seasonal dynamics, likely reflecting distinct local agricultural practices. Under heavy cloud contamination, NDVI values across all regions collapse toward zero, obscuring phenological signals and preventing reliable stage differentiation. These findings underscore the need for robust reconstruction techniques to restore NDVI time series integrity in cloud-prone regions.

### 2.4. A two-stage spatiotemporal reconstruction approach

Let the observation sites be indexed by $i(i = 1, ..., N)$ and the time steps by $t(t = 1, ..., T)$. For each site $i$ and time $t$, the input consists of three elements: the NDVI value $NDVI_{it}$, the CloudScore+ value $CS_{it}$, and the geographic coordinates $(LAT_i, LON_i)$. The proposed framework reconstructs cloud-contaminated NDVI time series through a two-stage process, spatiotemporal filling and multistep smoothing (Fig. 5).

#### 2.4.1. Stage 1 - spatiotemporal filling

Consecutive cloud-affected observations are reconstructed by pattern matching that leverages both the target site's historical temporal information and spatiotemporal patterns from other sites. This stage consists of four sequential steps:

1. Initialization of reconstruction parameters

<figure><img src="image_5.png" /><figcaption>Fig. 2. Rice growth phases observed through monthly ASF field surveys. Vegetative 2 corresponds to peak canopy greenness, whereas Land preparation reflects minimal vegetation cover (bare or recently ploughed fields) prior to the subsequent planting cycle.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Monthly percentage of NDVI observations severely affected by cloud and shadow interference (CloudScore+ ≤ 0.1) across Garut, Indramayu, and Pontianak from January 2021 to December 2022. Cloud interference exceeded 40% in most months, with a pronounced seasonal peak during the wet season, reaching its maximum of 88% in February 2022.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. NDVI distributions under less cloudy (CloudScore+ > 0.1, light blue) and cloud-contaminated (CloudScore+ ≤ 0.1, red) conditions across rice growth phases in (a) Garut, (b) Indramayu, and (c) Pontianak during 2021–2022. In Garut and Indramayu, less cloudy NDVI observations follow clear phenological cycles—rising during Vegetative 2 and declining through Harvest and Land Preparation. By contrast, the seasonal pattern in Pontianak is less distinct, likely reflecting different local practices. Under cloud-contaminated conditions, NDVI values remain near zero in all regions.</figcaption></figure>

a) CloudScore+ Threshold ($\tau_{CS}$): An NDVI value $NDVI_{it}$ is marked as a gap, denoted $NDVI_{it}^G$, if its corresponding CloudScore+ value satisfies $CS_{it} \le \tau_{CS}$.

b) Minimum Consecutive Gaps ($\tau_G$): Define

$$g_{it} = \begin{cases} 1, & \text{if } NDVI_{it}^G \\ 0, & \text{otherwise} \end{cases} \qquad (2)$$

For any contiguous subsequence ($t_s, t_s+1, ..., t_e$) with $t_e \le T$, the gap length for site i is $G_i = \sum_{t=t_e}^{t_s} g_{it}$. If $G_i \ge \tau_G$, the subsequence is classified as a large gap.

a) Adjacent Observations (A): The number of valid observations surrounding a large gap that are included to form the target interval.

2. Detection of large gaps and target interval formation

The algorithm scans the NDVI time series of all observation sites to identify large gaps based on previously defined reconstruction parameters. Once detected, each large gap is expanded into a target interval by including a specified number of adjacent valid observations. These adjacent points may lie on the preceding side, the following side, or both, depending on the availability of valid observations. This design enables the method to handle ongoing cropping cycles where post-gap data are not yet available, ensuring that reconstruction can still proceed using only preceding observations.

The length of a target interval is determined by the gap size ($G_i$) and the number of adjacent observations (A). If valid observations are available on both sides of the gap, the interval length equals $G_i + 2A$. If valid observations are available only on one side (either preceding or following), the interval length equals $G_i + A$ (Fig. 6a).

3. Identification of similar intervals

For each target interval, the algorithm performs a sliding-window
<figure><img src="image_8.png" /><figcaption>Fig. 5. Two-stage spatiotemporal reconstruction framework. Stage 1-Spatiotemporal Filling addresses large gaps by detecting target intervals and identifying matches using both the target site's historical temporal patterns and spatiotemporal patterns from other observation sites. Stage 2-Multistep Smoothing refines the reconstructed series through initial and iterative smoothing, aligning cloud-contaminated values with clearer observations.</figcaption></figure>

search across both the target site's historical temporal pattern and spatiotemporal patterns of other observation sites to identify similar temporal segments (see Fig. 6b). The degree of similarity is quantified using the Pearson correlation coefficient (Eq. 3):

$$ \rho = \frac{\sum (\mathbf{X} - \bar{\mathbf{X}})(\mathbf{Y} - \bar{\mathbf{Y}})}{\sqrt{\sum (\mathbf{X} - \bar{\mathbf{X}})^2 \sum (\mathbf{Y} - \bar{\mathbf{Y}})^2}} \quad (3) $$

Where X and Y represent the valid NDVI values from the target and candidate intervals, respectively. The candidate interval with the highest correlation, provided it exceeds 0.9, is retained for each observation site.

Candidate intervals were drawn from all available observation sites rather than being restricted to a fixed spatial radius. Each observation site corresponds to a 100 m × 100 m grid cell nested within a larger 300 m × 300 m block. While the ASF survey collected ground observations for the center pixel of every 100 m × 100 m grid cell, in this study, only one representative grid per 300 m × 300 m block was randomly selected for candidate searching. This choice was motivated by three considerations: (i) cloud conditions within a block are generally similar, reducing redundancy; (ii) random selection avoids systematic bias while ensuring diversity in candidate intervals; and (iii) limiting the search to one representative pixel per block improves computational efficiency without sacrificing reconstruction accuracy. In addition, the Euclidean distance between each candidate and the target site was calculated and incorporated into the reconstruction weighting scheme.

### 4. Candidate interval-based gap reconstruction

Each candidate interval is used to estimate the cloud-contaminated NDVI values within the large gap of the target interval (Fig. 6c)

through a third-degree polynomial regression (Eq. 4):

$$ \hat{\mathbf{Y}}_k = \beta_0 + \beta_1 \mathbf{X}_k + \beta_2 (\mathbf{X}_k)^2 + \beta_3 (\mathbf{X}_k)^3 \quad (4) $$

Where $\hat{\mathbf{Y}}_k$ is the NDVI predicted from candidate $k$ ($k = 1, ..., K$), $X_k$ is the corresponding NDVI value, and $\beta_p$ ($p = 0, 1, 2, 3$) are the regression coefficients estimated via least squares.

The final reconstructed NDVI, $\hat{\mathbf{Y}}$, is obtained by aggregating predictions from all K candidates using a weighted average (Eq. 5):

$$ \hat{\mathbf{Y}} = \frac{\sum_{k=1}^{K} \rho_k \tilde{\mathbf{d}}_k \hat{\mathbf{Y}}_k}{\sum_{k=1}^{K} \rho_k \tilde{\mathbf{d}}_k} \quad (5) $$

Where $\rho_k$ is the Pearson correlation and $\tilde{\mathbf{d}}_k$ is the normalized inverse Euclidean distance (Eq. 6) between the target interval and candidate interval k:

$$ \tilde{\mathbf{d}}_k = \frac{d_k - \min(\mathbf{d})}{\max(\mathbf{d}) - \min(\mathbf{d})} \quad (6) $$

Here, $d_k$ denotes the inverse Euclidean distance between the target interval and candidate interval k, while $\mathbf{d}$ represents the set of inverse Euclidean distances between the target site and all candidate sites. This weighting scheme integrates both temporal similarity ($\rho_k$) and spatial proximity ($\tilde{\mathbf{d}}_k$), assigning greater weight to closer sites while still allowing contributions from more distant sites when nearby sites are also cloud-contaminated. Closer sites are generally more reliable references, as they often share similar climate, environmental conditions, and cropping cycles [2]. If no candidate interval exceeds the correlation threshold, reconstruction defaults to Stage 2 (Multistep Smoothing).

### 2.4.2. Stage 2 - multistep smoothing

This stage refines the output of Stage 1 by correcting residual noise and minor gaps, ensuring that the reconstructed series conforms to expected vegetation dynamics. It is based on the assumption that NDVI trajectories evolve smoothly over time, reflecting gradual changes in vegetation phenology [9,20]. The objective is to progressively replace cloud-contaminated values while retaining high-confidence observations as temporal anchors.

Cubic splines are employed due to their flexibility and robustness in curve fitting [17,18,37]. Unlike methods such as the Savitzky-Golay filter, cubic splines avoid manually tuning parameters such as window size and polynomial degree [26]. The smoothing function is expressed as (Eq. 7):

$$ \tilde{Y}_{it} = f(t) \quad (7) $$

Where $\tilde{Y}_{it}$ is the smoothed NDVI value at site i and time t, and $f(t)$ is a smooth function of time. The spline is estimated by minimizing the weighted penalized sum of squares (Eq. 8):

$$ \min_{f} \sum_{i=1}^{N} \sum_{t=1}^{T} W_{it} (Y_{it} - \tilde{Y}_{it})^2 + \lambda \int (f''(t))^2 dt \quad (8) $$

Where $W_{it}$ is the weight at site i and time t, $\lambda$ is the smoothing parameter, and $f''(t)$ penalizes excessive curvature. The optimal $\lambda$ is determined by minimizing the Generalized Cross-Validation (GCV) score [38], as shown in Eq. (9):

$$ V(\lambda) = \frac{NT \sum_{i=1}^{N} \sum_{t=1}^{T} (Y_{it} - \tilde{Y}_{it})^2}{(NT - tr(\mathbf{A}))^2} \quad (9) $$

Where $V(\lambda)$ is the GCV score and $\mathbf{A}$ is the influence (hat) matrix. This stage proceeds in two sequential steps:

#### 1. Initial Smoothing

A weighted cubic spline is fitted to the Stage 1 NDVI series using CloudScore+ values as weights ($W_{it}$). Observations with CloudScore+
<figure><img src="image_9.png" /><figcaption>Fig. 6. Illustration of the Spatiotemporal Filling algorithm for reconstructing large gaps in NDVI time series. (a) Large gaps are detected and expanded into target intervals, including cloud-contaminated NDVI values (red) and their adjacent valid observations (black). (b) Candidate intervals are identified from both the target site's historical temporal patterns and the spatiotemporal patterns from other observation sites (green). Candidate intervals are generated by sliding the target interval across time and computing the Pearson correlation (ρ); Intervals with high similarity (ρ > 0.9) are retained. (c) Cloud-contaminated values within the large gap are reconstructed (blue) by aggregating predictions from all retained candidates using a weighted averaging scheme that integrates both temporal similarity and spatial proximity.</figcaption></figure>

below an initial threshold (default = 0.1) are replaced with spline predictions. Weights are then updated based on the normalized absolute difference between consecutive iterations (Eq. 10):

$$W_{it} = \frac{D_{it} - \min(D_i)}{\max(D_i) - \min(D_i)} \qquad (10)$$

Where $D_{it} = |NDVI_{it}^{(r+1)} - NDVI_{it}^{(r)}|$ is the difference between the smoothed NDVI at iteration r + 1 and r for site i and time t, and $D_i$ is the set of all such differences at site i.

## 2. Iterative Smoothing

The CloudScore+ threshold is progressively increased (e.g., 0.1, 0.2, 0.3, ...). At each step, values below the threshold are replaced with spline predictions, and weights are updated using Eq. (10). This iterative refinement continues until the maximum threshold is reached, producing a final smoothed NDVI curve that preserves high-quality observations while enforcing temporal continuity.

## 2.5. Benchmark data reconstruction methods

Two widely used reconstruction methods were implemented as benchmarks to evaluate the proposed Two-Stage Spatiotemporal

Approach: the Savitzky-Golay filter and cubic spline smoothing. Frequency-domain approaches, such as Fourier-based models, were excluded due to their limitations in fitting asymmetric and irregular vegetation time series [8,19].

### 2.5.1. Savitzky-Golay filter

The Savitzky-Golay filter [34] is a standard method for smoothing vegetation index time series [39,42]. It operates by fitting low-degree polynomials within a moving window, effectively performing a weighted local regression [8]. The smoothed value $Y_j^*$ at position j is expressed as in Eq. (11):

$$Y_j^* = \frac{\sum_{l=-m}^{m} C_l Y_{j+l}}{N} \quad j = 1, \dots, n \qquad (11)$$

Where Y are the original values, $C_l$ are smoothing coefficients, and $N = 2m + 1$ is the window width, with m representing half the window span. A larger window results in stronger smoothing but may obscure short-term variations. In this study, a window size of 11 points was selected to balance smoothness and temporal detail.

### 2.5.2. Cubic spline smoothing

Cubic spline smoothing is a flexible interpolation method that fits piecewise cubic polynomials between observed points, ensuring
continuous first- and second-order derivatives at the knots [16]. Unlike high-order global polynomial regression, cubic splines adapt to local temporal dynamics and reduce overfitting risk. The method is controlled by a smoothing parameter $\lambda$ [6], which balances fit and smoothness by minimizing the penalized sum of squares as in Eq. (8):

## 2.6. Evaluation and validation of reconstruction methods

Accurate NDVI reconstruction is critical for reliable crop monitoring and growth-phase analysis. Both controlled simulations and real-world validations were conducted to assess the effectiveness of the proposed Two-Stage Spatiotemporal Approach relative to benchmark methods. Simulations provided controlled conditions for method comparison, while ground-truth data from the ASF survey enabled validation under operational contexts.

### 2.6.1. Simulated data generation

Synthetic NDVI series were generated using a sine-based model, as shown in Eq. (12):

$$f(x) = A \times \sin\left(\frac{2\pi}{B}(x - C)\right) + D \quad (12)$$

Where A denotes the amplitude, B the period, C the phase shift, and D the vertical offset. Parameters were set as follows: A = 0.3, B = 100, C ~ U(-50, 50), and D ~ U(-0.05, 0.05). Each synthetic series spanned 365 days at 50 sites, with longitude drawn from U(100, 102) and latitude from U(-3, -1). Observations were sampled every five days, and initial CloudScore+ values were drawn from U(0.9, 0.95) to simulate clear-sky conditions.

### 2.6.2. Evaluation scenarios

Two cloud-contamination scenarios were designed:

1. Random Cloud Cover: cloudy days were selected randomly at three levels, low (10%), moderate (30%), and high (50%).

2. Consecutive Cloud Cover: the same contamination levels were imposed on consecutive days to mimic prolonged cloudy periods.

In both scenarios, cloud-contaminated NDVI values were replaced with random draws from U(0.01, 0.1), and CloudScore+ values from U(0.01, 0.1), simulating heavy cloud interference.

### 2.6.3. Validation metrics for simulations

Reconstruction accuracy was quantified using two standard indicators: Root Mean Squared Error (RMSE; Eq. 13) and Mean Absolute Error (MAE; Eq. 14).

$$RMSE = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (Y_i - \hat{Y}_i)^2} \quad (13)$$

$$MAE = \frac{1}{N} \sum_{i=1}^{N} |Y_i - \hat{Y}_i| \quad (14)$$

Where $Y_i$ and $\hat{Y}_i$ denote the original cloud-free and reconstructed NDVI values, respectively. Both metrics were calculated by comparing reconstructed series against the original cloud-free simulated series, thereby enabling a systematic evaluation of each method under varying cloud conditions.

### 2.6.4. Ground-truth validation with rice crop phases

The reconstructed NDVI series was validated against rice growth phase observations from the ASF survey to assess operational relevance. Four datasets were compared: (1) original NDVI, (2) NDVI reconstructed using the proposed Two-Stage Spatiotemporal Approach, (3) Savitzky-Golay, and (4) cubic spline smoothing.

A random forest classifier [5] was employed due to its robustness and ability to capture non-linear vegetation dynamics. The same predictive features were used across all datasets: current NDVI values, temporal gradients, and lagged NDVI values, which jointly represent crop status and short-term evolution.

The dataset was randomly partitioned into training (80%) and testing (20%) subsets, and the classification was repeated 100 times to reduce variability from random partitioning. Performance was assessed using precision (Eq. 15), recall (Eq. 16), and F1 score (Eq. 17):

$$Precision = \frac{TP}{TP + FP} \quad (15)$$

$$Recall = \frac{TP}{TP + FN} \quad (16)$$

$$F_1 = \frac{2 \cdot (Precision \times Recall)}{Precision + Recall} \quad (17)$$

where TP, FP, and FN denote true positives, false positives, and false negatives, respectively. These metrics comprehensively evaluate how well reconstructed NDVI time series support accurate rice growth stage classification.

## 3. Results

### 3.1. Performance on simulated data

The reconstruction performance of the Two-Stage Spatiotemporal Method was compared against two benchmark techniques, Savitzky-Golay filtering and cubic spline smoothing, using simulated NDVI datasets under two cloud-contamination scenarios: randomly distributed cloudy days and consecutively clustered cloudy days (Fig. 7).

#### 3.1.1. Randomly distributed cloudy days

As shown in Fig. 7a–7f, both Savitzky-Golay and cubic spline smoothing exhibited a sharp decline in accuracy as cloud contamination increased. At 10% random cloud cover, RMSE and MAE values were relatively low, but errors rose substantially at 30% and 50%, highlighting their sensitivity to missing observations. In contrast, the Two-Stage Spatiotemporal Method maintained consistently low RMSE and MAE across all levels of random contamination, indicating resilience against scattered cloud interference.

#### 3.1.2. Consecutive cloudy days

Results under consecutive cloud cover (Fig. 7g–7l) show that at 10% contamination, all methods performed similarly. However, as cloudy periods extended to 30% and 50% of the series, errors from the Savitzky-Golay and cubic spline methods increased sharply. In contrast, the Two-Stage Spatiotemporal Method maintained substantially lower RMSE and MAE values, demonstrating strong robustness even during prolonged cloud-induced data gaps.

### 3.2. Performance on real NDVI data

After validation on simulated datasets, the Two-Stage Spatiotemporal Method was applied to real-world NDVI observations. A representative case study was selected from the Garut Region during the 2021–2022 rice growing season. This site was characterized by severe cloud contamination and distinct seasonal NDVI dynamics, making it an ideal setting to evaluate reconstruction performance under practical conditions.

#### Stage 1. Spatiotemporal filling

Reconstruction Parameters were selected through preliminary analysis to optimize the filling process:
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Reconstruction performance of three methods: Savitzky-Golay (green), Cubic Spline Smoothing (purple), and the Two-Stage Spatiotemporal Method (red), evaluated using RMSE and MAE on simulated NDVI data. Results are shown under two cloud contamination scenarios: (a–f) random cloudy days and (g–l) consecutive cloudy days.</figcaption></figure>

* CloudScore+ Threshold: 0.1, ensuring that NDVI values above this cutoff retained enough information to distinguish rice growth phases.

* Minimum Consecutive Gaps: 7, targeting extended cloudy periods.

* Adjacent Observations: 24, corresponding to approximately four months of Sentinel-2 acquisitions at 5-day intervals, which reflects the average length of a single rice cropping cycle in Indonesia [2].

Under these settings, the site exhibited gaps of up to 11 consecutive cloudy observations. NDVI values with CloudScore+ ≤ 0.1 (light blue, Fig. 8a) dropped to near-zero, reflecting severe cloud/shadow contamination. After applying spatiotemporal filling, large gaps (red, Fig. 8b) were reconstructed using temporally similar patterns from historical seasons and spatiotemporal patterns from other observation sites, while unaffected observations (dark green) were preserved.

**Stage 2. Multistep smoothing**

The filled series was further refined with a weighted natural cubic spline. In the first iteration (threshold = 0.1), only low-quality observations were replaced with spline predictions (red, Fig. 9a), while high-quality data (dark green) were retained. In subsequent iterations (Fig. 9b–9j), the threshold was incremented by 0.1 at each step (up to 1.0), progressively correcting increasingly contaminated observations while leaving more reliable points untouched.

This iterative refinement effectively realigned cloud-affected values with clearer observations, producing a smooth and coherent NDVI trajectory across the season. By anchoring the process to high-quality data and progressively incorporating noisier points, the method eliminated cloud-related distortions while preserving the subtle temporal dynamics of rice phenology.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Application of the Spatiotemporal Filling stage to an NDVI time series from the Garut region during the 2021–2022 season. (a) Original NDVI observations, where values under cloudy conditions (CloudScore+ ≤ 0.1, light blue) show near-zero values, while clearer observations (CloudScore+ > 0.1, dark green) capture rice growth dynamics. (b) The NDVI series was reconstructed after applying spatiotemporal filling, with large gaps (red) replaced using temporally similar patterns from historical seasons and other observation sites, while unaffected observations remain unchanged.</figcaption></figure>

### 3.2.1. Application to all observation sites

The Two-Stage Spatiotemporal Method was systematically applied to 2,460 sites across Garut, Indramayu, and Pontianak (Fig. 9). Implementation was conducted in Python within a Jupyter Notebook environment and executed on standard hardware (Intel Core i7 processor, 64-bit architecture, 16 GB RAM), highlighting its feasibility for operational use. The source code is available at the GitHub repository: https://github.com/byususen/stms.

During Stage 1 (Spatiotemporal Filling), 1,446 sites contained large gaps, of which 1,220 (≈84%) were successfully reconstructed in ~7 hours. Stage 2 (Multistep Smoothing) required ~25 min. per iteration, totaling ~4 hours. Overall, the complete workflow was finalized within ~11 hours, demonstrating scalability and computational efficiency for large-area applications.

The reconstructed NDVI profiles (Fig. 10) showed substantial improvements compared to the original cloud-contaminated series (Fig. 4). In Garut (Fig. 10a) and Indramayu (Fig. 10b), reconstructed NDVI curves captured rice phenology accurately: NDVI rose during Vegetative 1, peaked at Vegetative 2, and declined through the Generative, Harvest, and Land Preparation phases. In Pontianak (Fig. 10c), profiles remained relatively high even during harvest and land preparation, reflecting

continuous vegetation cover and multi-harvest practices typical of the region.

### 3.2.2. Validation through prediction performance

To evaluate the practical effectiveness of reconstruction, a random forest classifier was applied to four datasets: (1) Original NDVI series (light blue), (2) NDVI reconstructed with the Two-Stage Spatiotemporal Method (red), (3) Cubic spline smoothing (purple), and (4) Savitzky-Golay filter (green). Classification performance was assessed across all rice growth stages (Fig. 11; Table 1).

### 3.2.3. Precision

The Two-Stage Spatiotemporal Method achieved the highest precision in every growth phase, with improvements of up to +19% over the original NDVI. Gains were most pronounced in the Generative phase (Fig. 11a). By contrast, cubic spline smoothing improved Vegetative 1 and Land Preparation but reduced precision in Vegetative 2 and Generative, while the Savitzky-Golay filter showed only marginal improvement in Land Preparation and declines in later phases.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Multistep Smoothing applied to NDVI observations after the Spatiotemporal Filling stage. (a) Initial smoothing with a CloudScore+ threshold of 0.1, where only heavily cloud-contaminated values (grey) are replaced by spline predictions (red), while high-quality observations (green) are preserved. (b–j) Subsequent iterations progressively increase the threshold (0.2 to 1.0), correcting additional contaminated values at each step while retaining reliable observations. This iterative refinement produces a continuous NDVI trajectory that captures rice growth dynamics across the season.</figcaption></figure>

### 3.2.4. Recall

The two-stage method also delivered the strongest recall, improving by up to +24% compared to the original series (Fig. 11b). Cubic spline smoothing provided modest recall gains in Vegetative 1 and Land Preparation, while Savitzky-Golay contributed minor benefits in Land Preparation. However, both benchmarks performed worse than the Two-Stage Method in later phases.

### 3.2.5. F1-Score and overall accuracy

As summarized in Table 1, the Two-Stage Method provided the most consistent F1-score improvements, with gains of up to +21%. The

highest F1-scores were recorded in the Generative and Harvest phases, reflecting a strong balance between precision and recall. Substantial improvement was also seen in Land Preparation, where baseline precision was low, enhancing the ability to detect early-season transitions. Regarding overall classification accuracy, the Two-Stage Method reached 0.673 (+9% compared to the original NDVI at 0.618). In contrast, cubic spline smoothing yielded only marginal gains (0.620), while Savitzky-Golay slightly reduced accuracy.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Reconstructed NDVI values across rice growth stages following application of the Two-Stage Spatiotemporal Method during the 2021–2022 season. Boxplots compare less cloudy observations (CloudScore+ > 0.1; light blue) with cloud-contaminated observations (CloudScore+ ≤ 0.1; red) for three study regions: (a) Garut, (b) Indramayu, and (c) Pontianak.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Precision (a) and Recall (b) of rice growth phase classification using a random forest model. Results are shown for the original NDVI dataset and three reconstruction methods: Savitzky-Golay, cubic spline smoothing, and the Two-Stage Spatiotemporal Method.</figcaption></figure>

# 4. Discussion

## 4.1. Efficiency of the interval-based strategy

A key innovation of the framework is its interval-based reconstruction. Unlike methods that depend on full-time series or multi-year datasets, the proposed approach operates on short, localized intervals that capture the most relevant temporal dynamics. This strategy enhances computational efficiency and improves adaptability to heterogeneous planting schedules and diverse field practices common in

tropical regions such as Indonesia.

Nevertheless, certain transitional phases, most notably Land Preparation, remain difficult to classify due to spectral ambiguity. Low NDVI values during this phase can resemble bare soil, flooding, or residue cover, overlapping with both preceding and subsequent growth stages. Local practices such as leaving crop residues further blur spectral distinctions. These challenges highlight the need for complementary indicators or ancillary datasets to refine phase delineation.
<caption>**Table 1**</caption>

F1-scores and overall accuracy of rice growth phase classification using a random forest model. Results are shown for the original NDVI dataset and reconstructed series with Savitzky-Golay, cubic spline smoothing, and the proposed Two-Stage Spatiotemporal Method.

<table><thead><tr><th>Growth Phase/ Overall Accuracy</th><th>Original Data</th><th>Savitzky-Golay</th><th>Cubic Spline Smoothing</th><th>Two-stage Spatiotemporal</th></tr></thead><tbody><tr><td><strong>F1 Score</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Vegetative 1</td><td>0.508</td><td>0.501</td><td>0.564*</td><td><strong>0.580*</strong></td></tr><tr><td>Vegetative 2</td><td>0.556</td><td>0.550</td><td>0.538</td><td><strong>0.656*</strong></td></tr><tr><td>Generative</td><td>0.699</td><td>0.677</td><td>0.630</td><td><strong>0.769*</strong></td></tr><tr><td>Harvest</td><td>0.708</td><td>0.698</td><td>0.703</td><td><strong>0.714*</strong></td></tr><tr><td>Land Preparation</td><td>0.460</td><td>0.493*</td><td>0.537*</td><td><strong>0.558*</strong></td></tr><tr><td><strong>Accuracy</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Overall</td><td>0.618</td><td>0.613</td><td>0.620</td><td><strong>0.673*</strong></td></tr></tbody></table>

Note: *Statistically significant improvements (95% level) over the original data

## 4.2. Predictive accuracy and robustness

The proposed method consistently outperformed benchmark techniques in both simulated and real-world settings, yielding higher reconstruction accuracy and improved rice growth phase classification. Its robustness stems from two complementary mechanisms: (i) interval-based filling that leverages historical temporal segments and spatiotemporal patterns from other sites, and (ii) multistep smoothing that progressively refines reconstruction under variable cloud conditions.

However, the method depends on the availability of correlated reference sites, and performance declines when no suitable patterns exist. Polynomial regression was selected for its simplicity and suitability for short intervals with limited data, where only a small number of adjacent observations are available. Its computational efficiency makes it practical for large-scale applications. Nonetheless, polynomial regression can occasionally produce out-of-range NDVI values that require threshold correction. These trade-offs highlight both the practicality of the current design and opportunities for methodological refinement in future work.

## 4.3. Transferability across indices, crops, and climates

The reconstructed NDVI profiles demonstrated strong alignment with rice phenology in Garut and Indramayu, where seasonal transitions between growth phases were clearly captured. In contrast, profiles in Pontianak remained consistently elevated, even during harvest and land preparation, reflecting continuous vegetation cover and multi-harvest practices typical of equatorial systems. These results highlight the method's adaptability across diverse cropping systems and agro-climatic conditions.

Although demonstrated with NDVI, the framework is inherently transferable to other vegetation indices such as the Enhanced Vegetation Index (EVI) and the Soil-Adjusted Vegetation Index (SAVI) [27]. Its modular design, spatiotemporal pattern matching combined with iterative smoothing, can also be tailored to crops with different growth cycles. Moreover, while developed to address persistent cloud cover in the humid tropics, the approach is equally relevant in temperate regions where cloud interference is seasonal or intermittent, broadening its utility for agricultural monitoring globally.

## 4.4. Computational and operational considerations

The proposed method provides an efficient reconstruction solution that is well-suited for operational monitoring at regional scales. However, scaling to national or continental levels introduces runtime and data volume challenges. Parallelization and GPU acceleration represent viable strategies to address these computational bottlenecks, but their integration requires further development. For large-scale adoption,

streamlined workflows that combine computational enhancements with automated data pipelines will be critical to ensure both efficiency and practicality in operational settings.

## 4.5. Broader implications for agricultural monitoring

Enhanced NDVI reconstruction directly strengthens agricultural monitoring, especially in regions where persistent cloud cover reduces the reliability of satellite observations. The proposed method improves satellite data's temporal continuity, increases crop growth monitoring precision, enhances yield forecasting, and supports more informed agricultural policy and food security planning.

It should be emphasized that this study represents a preliminary step toward operational monitoring. In Indonesia, where rice plots are predominantly managed under successive “rice-on-rice” cycles with minimal crop rotation, the key challenge lies in reliably identifying rice plots. Once these plots are detected, the proposed framework effectively reconstructs NDVI trajectories, even in regions frequently affected by cloud cover. Importantly, because the algorithm can operate using only preceding clear-sky observations, it can be applied without waiting for subsequent data. This design makes the method adaptable to near-real-time implementation, offering a practical pathway toward in-season monitoring in cloud-prone tropical environments.

## 4.6. Future work

Future improvements could further expand the utility and scalability of this method. Integrating additional satellite platforms and vegetation indices would extend its robustness across diverse remote sensing products. Application to a broader range of crops and agricultural systems would ensure adaptability beyond rice. Incorporating Synthetic Aperture Radar (SAR) data is a particularly promising direction, as radar penetrates cloud cover and can complement optical data in persistently cloudy regions. Finally, enhancing computational efficiency through parallelization or GPU-based processing will be essential for scaling to national and continental monitoring programs. Together, these extensions would broaden the framework's role in precision agriculture, enhance food system resilience, and enable operational large-scale crop monitoring.

# 5. Conclusion

This study presents a two-stage spatiotemporal reconstruction framework designed to address persistent cloud and shadow contamination in the NDVI time series, particularly relevant to tropical agricultural systems such as Indonesia. By integrating spatiotemporal pattern matching from historical seasons and other observation sites with a multistep temporal smoothing procedure, the method effectively reconstructs large observational gaps while adapting to heterogeneous planting schedules and regional agricultural practices.

Validation using simulated datasets showed that the proposed method consistently achieved the lowest RMSE and MAE compared to benchmark techniques, confirming its superior reconstruction performance. Application to real Sentinel-2 NDVI data for rice growth phase classification further demonstrated substantial gains in precision, recall, and F1-score relative to widely used alternatives such as the Savitzky-Golay filter and cubic spline smoothing.

These results demonstrate the proposed approach's robustness, scalability, and operational value. By reliably recovering cloud-contaminated observations and improving classification accuracy under cloud-prone conditions, the framework provides a practical solution for agricultural monitoring. Beyond methodological advancement, it contributes to global food security efforts by enabling more accurate yield assessments, supporting evidence-based policy planning, and strengthening applications in precision agriculture.
## CRediT authorship contribution statement

**B. Suseno:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**G. Brunel:** Writing – review & editing, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.
**H. Wijayanto:** Validation, Supervision, Methodology, Funding acquisition, Conceptualization.
**K. Sadik:** Supervision, Methodology, Conceptualization.
**F.M. Afendi:** Supervision, Methodology, Conceptualization.
**B. Tisseyre:** Writing – review & editing, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Bayu Suseno reports financial support and travel were provided by Erasmus. Bayu Suseno reports financial support was provided by IFI, French Embassy in Indonesia. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgment

This research was supported by l'Ambassade de France and l'Institut Français d'Indonésie under the Séjour Scientifique de Haut Niveau (SSHN) 2024 program (No. 165161Z); the Erasmus+ program through the European project MIC 2022 - KA171-77785; and scholarships awarded by BPS-Statistics Indonesia (Agreement No. 110/PTB-01/VIII/ 2021). We gratefully acknowledge their support, which made this study possible.

## Data availability

The authors do not have permission to share the data.

## References

[1] N. Agustiani, I. Gunawan, S. Margaret, S. Sujinah, Pola tanam padi untuk produktivitas tinggi dan indeks pertanaman yang optimal di lahan rawa pasang surut, Int. J. Geogr. Inf. Sci. 50 (3) (2022) 257–265, https://doi.org/10.24831/jai.v50i3.40983.

[2] M. Ardiansyah, H. Wijayanto, A. Kurnia, A. Djuraidah, Kajian perbandingan rata-rata umur panen dan produktivitas padi antar-provinsi di Indonesia, J. Penelit. Pertan. Tanam. Pangan 6 (1) (2022) 35, https://doi.org/10.21082/jpptp.v6n1.2022.p35-40.

[3] M. Ardiansyah, H. Wijayanto, A. Kurnia, A. Djuraidah, Multiclass forecasting on panel data using autoregressive multinomial logit and c5.0 decision tree, Pak. J. Stat. Oper. Res. (2023) 145–154, https://doi.org/10.18187/pjsor.v19i1.4053.

[4] P.M. Atkinson, C. Jeganathan, J. Dash, C. Atzberger, Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology, Remote Sens. Environ. 123 (2012) 400–417, https://doi.org/10.1016/j.rse.2012.04.001.

[5] L. Breiman, Random forests, Mach. Learn. 45 (1) (2001) 5–32, https://doi.org/10.1023/A:1010933404324.

[6] Z. Cai, P. Jönsson, H. Jin, L. Eklundh, Performance of smoothing methods for reconstructing ndvi time-series and estimating vegetation phenology from MODIS data, Remote Sens. (Basel) 9 (12) (2017) 12, https://doi.org/10.3390/rs9121271.

[7] R. Cao, Y. Chen, M. Shen, J. Chen, J. Zhou, C. Wang, W. Yang, A simple method to improve the quality of NDVI time-series data by integrating spatiotemporal information with the Savitzky-Golay filter, Remote Sens. Environ. 217 (2018) 244–257, https://doi.org/10.1016/j.rse.2018.08.022.

[8] J. Chen, Per. Jönsson, M. Tamura, Z. Gu, B. Matsushita, L. Eklundh, A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter, Remote Sens. Environ. 91 (3) (2004) 332–344, https://doi.org/10.1016/j.rse.2004.03.014.

[9] D. Chu, H. Shen, X. Guan, J.M. Chen, X. Li, J. Li, L. Zhang, Long time-series NDVI reconstruction in cloud-prone regions via spatio-temporal tensor completion,

Remote Sens. Environ. 264 (2021) 112632, https://doi.org/10.1016/j.rse.2021.112632.

[10] Delince, J., Lemoine, G., Defourny, P., Gallego, J., Davidson, A., Shibendu Ray, Rojas, O., Latham, J., & Achard, F. (2017). Handbook on remote sensing for agricultural statistics. 10.13140/RG.2.2.132596.9920.

[11] M. Drusch, U. Del Bello, S. Carlier, O. Colin, V. Fernandez, F. Gascon, B. Hoersch, C. Isola, P. Laberinti, P. Martimort, A. Meygret, F. Spoto, O. Sy, F. Marchese, P. Bargellini, Sentinel-2: Esa's optical high-resolution mission for games operational services, Remote Sens. Environ. 120 (2012) 25–36, https://doi.org/10.1016/j.rse.2011.11.026.

[12] European Space Agency. (2022). Sentinel-2 msi level-2a boa reflectance [dataset]. https://doi.org/10.5270/S2-znk9xsj.

[13] L. Gandharum, M.E. Mulyani, D.M. Hartono, A. Karsidi, M. Ahmad, Remote sensing versus the area sampling frame method in paddy rice acreage estimation in Indramayu regency, west java province, Indonesia, Int. J. Remote Sens. 42 (5) (2021) 1738–1767, https://doi.org/10.1080/01431161.2020.1842541.

[14] F. Gao, X. Zhang, Mapping crop phenology in near real-time using satellite remote sensing: challenges and opportunities, J. Remote Sens. 2021 (2021) 8379391, https://doi.org/10.34133/2021/8379391, 2021/.

[15] S. Ghildiyal, N. Goel, S. Singh, S. Lal, R. Kawsar, A.E. Saddik, M. Saini, SSGAN: cloud removal in satellite images using spatiospectral generative adversarial network, Eur. J. Agron. 161 (2024) 127333, https://doi.org/10.1016/j.eja.2024.127333.

[16] P.J. Green, Bernard.W. Silverman, Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach, 0, Chapman and Hall/CRC, 1993, https://doi.org/10.1201/b15710.

[17] T.J. Hastie, R.J. Tibshirani, Generalized Additive Models, 1st ed., Routledge, 2017 https://doi.org/10.1201/9780203753781.

[18] T. Hastie, R. Tibshirani, Generalized additive models, Stat. Sci. 1 (3) (1986), https://doi.org/10.1214/ss/1177013604.

[19] J.N. Hird, G.J. McDermid, Noise reduction of NDVI time series: an empirical comparison of selected techniques, Remote Sens. Environ. 113 (1) (2009) 248–258, https://doi.org/10.1016/j.rse.2008.09.003.

[20] Y. Julien, J.A. Sobrino, Comparison of cloud-reconstruction methods for time series of composite NDVI data, Remote Sens. Environ. 114 (3) (2010) 618–625, https://doi.org/10.1016/j.rse.2009.11.001.

[21] Y. Julien, J.A. Sobrino, Optimizing and comparing gap-filling techniques using simulated NDVI time series from remotely sensed global data, Int. J. Appl. Earth Obs. Geoinf. 76 (2019) 93–111, https://doi.org/10.1016/j.jag.2018.11.008.

[22] I.N. Khasanah, D.W. Triscowati, A. Handoyo, W.P. Buana, K. Kadir, Kombinasi survei kerangka sampel area dan remote sensing untuk estimasi luas tanaman padi di masa pandemi (ksa-hybrid): studi kasus di provinsi lampung, J. Appl. Komputasi Stat. 12 (2) (2020) 53, https://doi.org/10.34123/jurnalasks.v12i2.284.

[23] Y. Kurniawati, H. Wijayanto, A. Kurnia, D. Dirgahayu D, B. Susetyo, Rice phenology monitoring via ensemble classification for an extremely imbalanced multiclass dataset of hybrid remote sensing, Remote Sens. Appl.: Soc. Environ. 35 (2024) 101246, https://doi.org/10.1016/j.rsase.2024.101246.

[24] D. Kusumaningrum, H. Wijayanto, A. Kurnia, K.A. Notodiputro, M. Ardiansyah, I. M. Parvez, Four-parameter beta mixed models with survey and sentinel 2A satellite data for predicting paddy productivity, Smart Agric. Technol. 9 (2024) 100525, https://doi.org/10.1016/j.atech.2024.100525.

[25] S. Li, L. Xu, Y. Jing, H. Yin, X. Li, X. Guan, High-quality vegetation index product generation: a review of NDVI time series reconstruction techniques, Int. J. Appl. Earth Obs. Geoinf. 105 (2021) 102640, https://doi.org/10.1016/j.jag.2021.102640.

[26] Y. Liu, B. Dang, Y. Li, H. Lin, H. Ma, Applications of savitzky-golay filter for seismic random noise reduction, Acta Geophys. 64 (1) (2016) 101–124, https://doi.org/10.1515/acgeo-2015-0062.

[27] L. Lu, C. Kuenzer, C. Wang, H. Guo, Q. Li, Evaluation of three modis-derived vegetation index time series for dryland vegetation dynamics monitoring, Remote Sens. (Basel) 7 (6) (2015) 6, https://doi.org/10.3390/rs70607597.

[28] R. Michishita, Z. Jin, J. Chen, B. Xu, Empirical comparison of noise reduction techniques for NDVI time-series based on a new measure, ISPRS J. Photogramm. Remote Sens. 91 (2014) 17–28, https://doi.org/10.1016/j.isprsjrs.2014.01.003.

[29] V.J. Pasquarella, C.F. Brown, W. Czerwinski, W.J. Rucklidge, Comprehensive quality assessment of optical satellite imagery using weakly supervised video learning, in: Proceedings of the 023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), 2023, pp. 2125–2135, https://doi.org/10.1109/CVPRW59228.2023.00206.

[30] F. Rodríguez-Puerta, R.L. Perroy, C. Barrera, J.P. Price, B. García-Pascual, Five-year evaluation of sentinel-2 cloud-free mosaic generation under varied cloud cover conditions in Hawai'i, Remote Sens. (Basel) 16 (24) (2024) 24, https://doi.org/10.3390/rs16244791.

[31] G.J. Roerink, M. Menenti, W. Verhoef, Reconstructing cloudfree NDVI composites using fourier analysis of time series, Int. J. Remote Sens. 21 (9) (2000) 1911–1917, https://doi.org/10.1080/014311600209814.

[32] J.W. Rouse, R.H. Haas, J.A. Schell, D.W. Deering, Monitoring vegetation systems in the great plains with ERTS. https://ntrs.nasa.gov/citations/19740022614, 1974.

[33] K. Ruslan, Improving Indonesia's Food Statistics Through The Area Sampling Frame Method, 0, Center for Indonesian Policy Studies, 2019, https://doi.org/10.35497/287781.

[34] Abraham. Savitzky, M.J.E. Golay, Smoothing and differentiation of data by simplified least squares procedures, Anal. Chem. 36 (8) (1964) 1627–1639, https://doi.org/10.1021/ac60214a047.
[35] Jean. Steinier, Yves. Termonia, Jules. Deltour, Smoothing and differentiation of data by simplified least square procedure, Anal. Chem. 44 (11) (1972) 1906–1909, https://doi.org/10.1021/ac60319a045.

[36] UN.ESCAP, Using big data for official statistics: key considerations when using mobile phone data. https://repository.unescap.org/handle/20.500.12870/5177, 2022.

[37] N. Vorobiova, A. Chernov, Curve fitting of modis ndvi time series in the task of early crops identification by satellite images, Procedia Eng. 201 (2017) 184–195, https://doi.org/10.1016/j.proeng.2017.09.596.

[38] S.N. Wood, Generalized Additive Models: An Introduction with R, 2nd ed., Chapman and Hall/CRC, 2017 https://doi.org/10.1201/9781315370279.

[39] F. Xie, H. Fan, Deriving drought indices from MODIS vegetation indices (NDVI/ EVI) and land surface temperature (LST): is data reconstruction necessary? Int. J.

Appl. Earth Obs. Geoinf.Int. J. Appl. Earth Obs. Geoinf. 101 (2021) 102352
https://doi.org/10.1016/j.jag.2021.102352.

[40] L. Zeng, B.D. Wardlow, D. Xiang, S. Hu, D. Li, A review of vegetation phenological metrics extraction using time-series, multispectral satellite data, Remote Sens. Environ. 237 (2020) 111511, https://doi.org/10.1016/j.rse.2019.111511.

[41] J. Zhang, J. Xiao, X. Tong, J. Zhang, P. Meng, J. Li, P. Liu, P. Yu, NIRv and SIF better estimate phenology than NDVI and EVI: effects of spring and autumn phenology on ecosystem production of planted forests, Agric. For. Meteorol. 315 (2022) 108819, https://doi.org/10.1016/j.agrformet.2022.108819.

[42] J. Zhou, L. Jia, M. Menenti, B. Gorte, On the performance of remote sensing time series reconstruction methods – a spatial comparison, Remote Sens. Environ. 187 (2016) 367–384, https://doi.org/10.1016/j.rse.2016.10.025.