A Multi-Stage Framework for Cloud-Free NDVI Generation in Agricultural Monitoring Using Sentinel-2 Imagery

Abstract

Cloud cover in optical satellite imagery is a major impediment to continuous agricultural monitoring, often disrupting the derivation of Normalized Difference Vegetation Index (NDVI) time series. To address this challenge, this study introduces a multi-stage framework to generate gap-free (cloud-free) NDVI images for agricultural monitoring using Sentinel-2 imagery. The proposed framework comprises three sequential modules: (1) Constrained Cloud-Maximum NDVI Composite (CC-MNC), which produces an initial NDVI composite with minimal cloud contamination; (2) local median-based spatial interpolation, which fills residual small-area gaps by replacing missing pixels with the median NDVI of nearby clear-sky pixels within a local neighborhood; and (3) temporal gap-filling using Gaussian Process Regression (GPR), which estimates NDVI values for larger missing regions or periods affected by persistent cloud cover. The framework was evaluated on five agricultural reference tiles in South Korea during a cloud-prone summer season. The initial Sentinel-2 images had an average cloud cover of 21.37%. After applying CC-MNC, the average cloud coverage in the NDVI composites was reduced to 0.93%. Subsequent spatial gap-filling further decreased the remaining missing data to 0.47%. Finally, the temporal GPR module eliminated all remaining gaps, achieving a fully cloud-free NDVI product (0% cloud cover) with a 100% gap recovery rate. The average processing time to generate each composite was under 10 minutes, demonstrating the method's efficiency for near-real-time applications. These results indicate that the proposed multi-stage approach can effectively provide complete and reliable NDVI maps under persistent cloud conditions. The method is readily applicable to future agricultural satellite missions,such as Compact Advanced Satellite 500-4 (CAS500-4), to enhance operational crop monitoring during cloudy seasons.

1. Introduction

Global warming and increasingly frequent extreme weather events have amplified the variability in crop growth timing, thereby increasing the importance of Earth observation satellite imagery that can continuously monitor crop conditions in space and time (Dalposso et al., 2013; Weiss et al., 2020; Li et al., 2021). In operational crop monitoring, the Normalized Difference Vegetation Index (NDVI) derived from optical satellite sensors is widely used, as it is a key vegetation index closely related to crop vigor, above-ground biomass, and yield. In particular, to capture rapid phenological transitions during the growing season—such as emergence, stem elongation, heading/flowering, and harvest—it is essential to secure a temporally dense and stable NDVI time series. However, optical sensors are inherently constrained by weather conditions, so minimizing the impact of clouds and other adverse atmospheric conditions is crucial for reliable crop monitoring.

Korea is characterized by frequent cloud occurrence due to the summer monsoon, typhoons, and spring–autumn haze and fog. These regional climatic conditions often result in missing satellite observations, posing a major obstacle to establishing a stable agricultural monitoring system (Roy et al., 2008; Moreno-Martínez et al., 2020). To mitigate this problem, a variety of NDVI compositing and gap-filling techniques, ranging from traditional and statistical approaches to machine-learning-based–based methods, have been proposed and applied worldwide (Gutman, 1999; Jung et al., 2012; Belda et al., 2020; Chu et al., 2021; Li et al., 2021; Garioud et al., 2021; Yu et al., 2021; Eun et al., 2023; Kim et al., 2024).

For example, Belda et al. (2020) applied Gaussian Process Regression (GPR) to Sentinel-2-based Leaf Area Index (LAI) time series and showed that using crop-type-specific hyperparameters can substantially improve computational efficiency without degrading accuracy. Chu et al. (2021) proposed an adaptive spatio-temporal tensor completion (ST-Tensor) method that simultaneously exploits spatial neighborhood, inter-annual variability, and seasonality to reconstruct long-term Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI series (2001 to 2018) over cloud-prone Southeast Asia, outperforming five existing gap-filling methods in both temporally continuous and spatio-temporally continuous gap scenarios. Yu et al. (2021) introduced a climate-integrated gap-filling (CGF) method that combines Landsat NDVI with ECMWF Re-Analysis (ERA) reanalysis climate data using a Radial Basis Function Network (RBFN) and demonstrated that it can reconstruct 30 m, 8 day NDVI time series with lower error than HANTS, Savitzky–Golay filtering, and climatology-based approaches, even in the presence of large temporal gaps.

However, studies on the generation and reconstruction of composite NDVI products that reflect the characteristics of Korean croplands are still limited. For example, the Maximum NDVI Composite (MaxNDVI) method is one of the most commonly used techniques for reducing the impact of clouds (Kim and Eun, 2022). MaxNDVI effectively suppresses cloud-contaminated observations within a given compositing window, but when the number of input images is limited, or clouds persist over an extended time, large gaps can still remain in the composite NDVI images (Griffiths et al., 2019). Conversely, purely temporal methods such as GPR can restore missing values along the time axis but are sensitive to residual cloud contamination or spatial artifacts in the input NDVI series and may propagate such errors to multiple composite periods. Such missing or contaminated areas can become a major source of error in downstream applications such as crop condition assessment, phenological stage estimation, and yield prediction (Belda et al., 2020; Eun et al., 2021; Kim et al., 2024).

In particular, Korean agricultural fields are typically small and spatially heterogeneous, which implies that an effective reconstruction strategy should combine spatial correction for small-scale gaps with temporal correction for longer cloud-induced gaps so as to preserve both field-level spatial patterns and temporal continuity (van Leeuwen et al., 1999; Wang et al., 2014; Li et al., 2021; Belda et al., 2020; Kim et al., 2024). The objective of this study is therefore to develop and present a gap-minimized NDVI production framework for Korean agricultural regions by integrating Constrained Cloud-Maximum NDVI Composite (CC-MNC), local median-based spatial gap-filling, and GPR-based temporal reconstruction based on Sentinel-2 imagery, and to quantitatively evaluate its performance. To this end, five standard grid tiles (Tile IDs: 96, 97, 109, 121, 130) located in agricultural areas were selected, and 15 day composite NDVI images were generated for each tile using three daily NDVI scenes and their corresponding cloud masks as inputs. Spatial and temporal gap-filling procedures were then sequentially applied to the composite images to assess how effectively cloud-induced missing areas could be reconstructed. In addition, considering operational system implementation, we defined metadata for the resulting products and quantitatively evaluated processing times under specific conditions.

2. Materials and Methods

2.1. Study Area

The C4SW Standard Grid is a tiling scheme defined in the Korea 2000 / Unified Coordinate System (Korea_2000_Korea_Unified_Coordinate_System, EPSG:5179), which is the national unified coordinate system officially adopted by the National Geographic Information Institute of Korea, and it serves as the standard grid system for the Korean agricultural satellite mission, specifying the basic tile units on which NDVI and other agricultural monitoring products will be generated (Kim et al., 2024). In this study, five tiles (Tile IDs: 96, 97, 109, 121, 130) were selected from the C4SW Standard Grid as representative units for major agricultural areas on the Korean Peninsula (Fig. 1). These tiles cover key paddy and upland cropping regions along the west coast and in inland areas, including rice paddies around Dangjin and upland fields in the central inland and southern regions, thereby capturing typical cropping environments in Korea.

Fig. 1. Study area (Tile IDs: 96, 97, 109, 121, 130; red hatched) selected within the C4SW Standard Grid over the Korean peninsula agricultural region.

2.2. Satellite Data

In this study, we additionally use NDVI composite products from MODIS, PlanetScope, and Landsat-8/9 as baseline datasets for comparison, in order to evaluate the performance of the CC-MNC composites generated from Sentinel-2 imagery.

2.2.1. MODIS

MODIS is a multi-spectral sensor aboard NASA’s Terra and Aqua satellites, designed for global Earth observation. It measures data in 36 spectral bands from the visible through shortwave infrared to thermal infrared, with spatial resolutions of approximately 250 m, 500 m, and 1,000 m depending on the band. With a wide viewing swath of about 2,330 km, MODIS can cover the entire globe every 1 to 2 days. The sensor’s data products include calibrated and atmospherically-corrected reflectance images and various derived indices, along with other parameters like land surface temperature and albedo. These characteristics make MODIS well-suited for long-term, broad-scale time-series analyses of vegetation and land surface dynamics.

2.2.2. PlanetScope

PlanetScope is a commercial constellation of numerous small satellites that provides near-daily optical imagery of almost all Earth’s land areas at approximately 3 to 4 m spatial resolution. Each PlanetScope satellite carries a sensor with four spectral bands with Blue, Green, Red, and Near-Infrared (NIR). Owing to the large number of satellites, the same location can be observed on an almost daily basis, and the data are delivered as orthorectified surface reflectance imagery with radiometric and atmospheric corrections applied. This high-frequency, high-resolution data is very useful for detailed spatial pattern analysis at the field or parcel scale, making PlanetScope imagery well-suited for applications like crop monitoring and fine-grained land cover mapping.

2.2.3. Landsat-8/9

Landsat-8/9 are moderate-resolution Earth observation satellites jointly operated by the United States Geological Survey (USGS) and NASA, continuing the decades-long Landsat program. Both satellites carry identical instruments: the Operational Land Imager (OLI) for visible, near-infrared, and shortwave infrared imaging, and the Thermal Infrared Sensor (TIRS) for thermal imaging. The OLI provides 30 m resolution multi-spectral imagery (with a 15 m panchromatic band), and TIRS acquires thermal infrared data at 100 m resolution which is resampled to 30 m in the standard products. Each satellite has a 16 day repeat orbital cycle; together, the two satellites offset each other to effectively revisit the same location approximately every 8 days. The data from Landsat-8/9 used in this study are provided as Level-2 Surface Reflectance products, which include radiometric and geometric corrections as well as atmospheric correction. These analysis-ready products are well-suited for long-term time-series studies and are widely used for land cover change detection and vegetation monitoring.

2.2.4. Sentinel-2A/B

We utilized imagery from the Sentinel-2A/B satellites operated by European Space Agency (ESA). In particular, Level-2A surface reflectance products with atmospheric correction applied were used to calculate NDVI. Sentinel-2 provides high spatial resolution with 10 m for the relevant red and near-infrared bands used in NDVI computation and a 5 day revisit cycle, making it well-suited for crop monitoring. For this study, we used Sentinel-2 images covering the study areas during the May through October, 2023 (Fig. 2).

Fig. 2. Sentinel-2 mosaic Red-Greed-Blue (RGB) image from tiles T52SBG, T52SCF, T52SBF, and T52SCG over the study region.

2.2.5. NDVI from Sentinel-2A/B Image

For this study, Sentinel-2 Level-2A surface reflectance products were used as the primary input for constructing the CC-MNC images. To mitigate the limitations of long compositing windows, the compositing period was restricted to 10 days, and a viewing angle constraint of Viewing Zenith Angle (VZA) ≤ 8.4° was applied. In the preprocessing step, each Sentinel-2 scene was (1) cloud-masked using the Scene Classification Layer (SCL), (2) clipped to the study area boundary defined by a shapefile, and (3) resampled to a common projection and grid system.

It should be noted that the Sentinel-2 SCL mask is not perfect. Thin clouds, cloud shadows, or bright surfaces such as snow and bare soil can occasionally be misclassified, leading to both omission and commission errors in the cloud mask. If such residual misclassifications are not handled properly, they may leave contaminated NDVI values that degrade the effectiveness of subsequent temporal gap-filling. For this reason, the proposed framework combines the SCL mask with additional spatial filtering and compositing steps to further suppress cloud-affected pixels before applying temporal reconstruction.

NDVI was then calculated from the 10 m bands using the Red (B04, 664 nm) and NIR (B08, 842 nm) bands as

\(\begin{align}N D V I =\frac{N I R - R e d}{N I R + R e d}\end{align}\)       (1)

2.3. Algorithms

The framework presented in this study comprises three algorithms. The initial input data is a daily vegetation index, and the process is carried out sequentially in the following order: periodic NDVI compositing, local median-based spatial gap-filling, and temporal gap-filling of missing data.

2.3.1. CC-MNC

First, to evaluate the periodic composite vegetation index generated from Sentinel-2 NDVI, we compared it with NDVI composite products derived from MODIS, Landsat-8, and PlanetScope imagery. NDVI for each sensor was calculated over two study areas, the Dangjin rice-paddy region and the Taebaek highland cabbage fields, for the period 2019–2021. Monthly maximum NDVI composites were then constructed and analyzed together with monthly precipitation to examine how differences in cloud cover and observation frequency affect the NDVI time series. In addition, NDVI derived from field spectroradiometer measurements collected in 2021 was used as a reference. The overall workflow of the multi-sensor comparison is summarized in Fig. 3, and the main findings from the NDVI and precipitation time-series analysis are described below (Eun et al., 2023).

Fig. 3. Flowchart of a study on the application of vegetation index-based CC-MNC image using various satellite data.

The 2019–2021 NDVI time-series comparison revealed several characteristic behaviors among the sensors. In 2020, when precipitation was exceptionally high, the Sentinel-2A/B NDVI series with a 5 day revisit showed large short-term fluctuations that reflect frequent cloud contamination. By contrast, PlanetScope NDVI, which is acquired only on relatively cloud-free dates, showed much weaker cloud-induced variability. The MODIS NDVI product, provided as a 16 day composite, effectively minimized the influence of clouds while still reproducing the seasonal growth pattern. Landsat-8 NDVI, on the other hand, remained systematically lower than the other satellite products throughout most of the season. Consistent with the spatial distribution, Sentinel-2A/B NDVI tended to be lower than the other sensors before reaching its seasonal maximum, and this tendency was particularly evident in the wet year 2020. PlanetScope and MODIS NDVI exhibited broadly similar temporal trajectories and were in good agreement with the field-derived NDVI. However, for MODIS, the coarse spatial resolution (250 m) caused NDVI to increase earlier and remain high for a longer period than for the other sensors, because reflectance from neighboring forests and upland vegetation was mixed into the MODIS pixels over paddy fields that still contained standing water in May (Eun et al., 2023).

Overall, MODIS NDVI reliably reproduced seasonal-scale vegetation growth variations, but its coarse spatial resolution meant that the influence of surrounding non-agricultural land cover was pronounced. PlanetScope NDVI, with meter-level spatial resolution, captured within-field heterogeneity and fine-scale spatial structures in detail; however, differences in radiometric characteristics among sensors, limited image acquisition frequency, and data costs constrain its use for constructing large-area, long-term NDVI time series. Landsat-8 NDVI offers moderate resolution (30 m) and the advantage of a long observation record, but its 16 day revisit interval can lead to data gaps in regions with persistent cloud cover. Although for the two example sites considered here, the Landsat-8 time series during the summer of 2020 appears more temporally continuous than Sentinel-2A/B because several clear-sky Landsat observations happened to be acquired on key dates, at the national scale, the two-satellite Sentinel-2 constellation provides 10 m spatial resolution and an approximately 5 day nominal revisit cycle, offering a higher potential sampling frequency under comparable cloud conditions. For this reason, Sentinel-2 NDVI was selected as the base dataset for generating periodic composite NDVI products over Korean agricultural areas. In addition, Sentinel-2 provides visible-to-NIR spectral bands, including red and near-infrared bands at 10 m resolution, that are broadly consistent with those planned for the Korean agricultural satellite mission, so the proposed compositing and gap-filling framework can be readily transferred to future operational NDVI products.

Because the viewing geometry of Sentinel-2 can vary for the same ground pixel, only observations with a viewing-incidence zenith angle less than or equal to 8.4°, as provided in the Level-2A metadata, were retained as composite candidates. For each 10 day period, the final composite NDVI image was generated using a MaxNDVI approach that selects, for each pixel, the maximum NDVI value among the eligible observations. In parallel, the residual cloud mask was not produced by simply aggregating daily SCL-based cloud masks within the compositing window. Instead, the SCL class of the observation selected by the MaxNDVI procedure was assigned to the corresponding composite pixel. Consequently, even if some observations for the same pixel were classified as cloud or cloud shadow, the pixel was not regarded as cloud in the composite result when a clear-sky observation existed within the window and was selected by MaxNDVI.

By exploiting the short revisit interval of the Sentinel-2A/B constellation, approximately 5 days, this procedure produced periodic NDVI composites with minimal missing data even during cloudy seasons (Figs. 4 and 5). The resulting 10 day Sentinel-2 MaxNDVI time series served as the input to the subsequent spatio-temporal gap-filling modules.

Fig. 4. CC-MNC period composite vegetation index production process prior to gap filling.

Fig. 5. Example of CC-MNC synthesis results using daily vegetation index data (Tile 110, 130).

2.3.2. Spatial Gap-Filling for Cloud Masking Boundaries

The second stage addresses any residual gaps in the composite images from Stage 1. Although CC-MNC greatly reduces cloud cover in the NDVI composites, there can still be small missing patches, especially in unusually cloudy months or in areas of complex terrain where cloud shadows or bright adjacent clouds might cause the algorithm to exclude some pixels. Additionally, if a shorter compositing period is used, more gaps are expected. Our gap-filling procedure comprises two sub-steps.

We first perform a spatial interpolation to fill in small isolated cloudy pixels or edges of cloud-masked regions within each composite image. For each pixel identified as missing in the CC-MNC process, we examine a surrounding window (e.g., 3 × 3 or 5 × 5 pixels) to find valid NDVI values of neighboring pixels. If at least half of the neighboring pixels are cloud-free, we assign the median of those neighbors’ NDVI values to the target pixel. This method assumes local spatial continuity of crop conditions and is effective for patching small holes from the boundary of a cloud mask. In cases where an entire area is missing, spatial filling is less reliable, so those are left for the temporal gap-filling step.

2.3.3. Temporal Gap-Filling Using GPR

For larger gaps that cannot be filled spatially—such as fields that are completely obscured by clouds in all input images—we apply temporal gap-filling using Gaussian Process Regression. For each pixel, an NDVI time series is constructed from the sequence of periodic composite images, and a Gaussian Process model is fitted to the available observations. The model then predicts the NDVI value for the missing composite period and simultaneously provides an estimate of the prediction uncertainty.

In this study, time is represented not in calendar days but as a discrete index of the composite periods of 10 day intervals, so that each step corresponds to one compositing window. We adopt a squared exponential kernel for the GPR, which is well-suited to representing the smooth seasonal trajectory of crop growth, and use the composite-period index as the input variable. When NDVI is missing at a target composite period T_K, the GPR model is trained using NDVI values from neighboring composite periods such as T_(K-2), T_(K-1), and T_(K+1), thereby interpolating or extrapolating the value at T_K based on both preceding and following observations within the same growing season (typically up to ±2 composite periods around the gap). While GPR provides a flexible and probabilistic framework for reconstructing missing NDVI values, its performance can degrade when observations are extremely sparse or unevenly distributed in time, and when outliers remain in the input series.

Moreover, pixel-wise GPR can be computationally expensive if applied directly to all daily observations over large areas. In the proposed system, these limitations are mitigated by applying GPR only to periodic composite NDVI after rigorous cloud masking and spatial gap-filling, thereby reducing noise in the input series and the number of time steps to be modeled. GPR is applied on a per-pixel basis for computational efficiency, and the resulting estimates preserve the spatial structure of fields while reconstructing only the cloud-contaminated areas. After temporal gap-filling, all pixels that were previously missing in the composite image are assigned physically plausible NDVI values, yielding a temporally continuous, gap-corrected NDVI field.

2.4. Framework

2.4.1. Workflow Overview

This section describes how the three modules introduced above are integrated into the operational processing chain of the Compact Advanced Satellite 500-4 (CAS500-4) satellite to produce the final NDVI products. The NDVI time-series reconstruction framework consists of four stages: (a) input data preparation, (b) CC-MNC-based periodic compositing, (c) spatial gap-filling, and (d) temporal gap-filling.

In the input data preparation stage, daily NDVI products (e.g., L3V NDVI) and the corresponding cloud masks provided by the satellite are collected and reorganized onto the standard observation grid (tile) system(Fig. 6a). Each daily image is mapped to the predefined tile coordinate system so that all input data are spatially consistent.

Fig. 6. Workflow of CC-MNC-based periodic NDVI compositing and two-stage spatio-temporal gap-filling: (a) data preparation, (b) CC-MNC periodic compositing, (c) spatial gap-filling, and (d) temporal gap-filling using GPR.

In the CC-MNC-based periodic compositing stage, daily NDVI and cloud mask information within each composite period are used to generate a representative NDVI for every tile and pixel (Fig. 6b). For each composite period, the CC-MNC algorithm selects an NDVI value using daily NDVI and cloud information accumulated over that period.

In the spatial gap-filling stage, a first spatial correction is applied to residual missing pixels that remain after CC-MNC compositing (Fig. 6c). The procedure targets small gaps and isolated cloudy pixels and fills them using spatially neighboring valid NDVI values.

Finally, in the temporal gap-filling stage, medium- to large-scale missing regions that still remain after spatial gap-filling are corrected using time-series information and a GPR model (Fig. 6d). For each pixel, a local NDVI time series is constructed from the composite NDVI at the target period together with the composites and cloud masks of the preceding and following periods. A GPR model is then fitted to this time series, and the NDVI value at the missing period is estimated from the fitted model. By inferring values that smoothly connect the NDVI changes before and after the gap, this step reconstructs the NDVI at the missing time. Through these stepwise corrections, the framework ultimately produces a tile-based periodic composite NDVI time series in which cloud-induced gaps have been removed.

2.4.2. Output Format and File Naming

The outputs generated by the proposed framework are defined as tile-based composite images in GeoTIFF format and an accompanying metadata file. The filename of each composite product follows the rule:

C4_{CompositeStartDate}_{CompositeEndDate}_{TileID}_L3VP_{ProductName}

A total of five products are provided: four image layers (NDVI, Normalized Difference Red Edge Index [NDRE], Flag, and Cloud) and one auxiliary metadata file (Meta(aux)). All image products are stored on a tile grid of 10,000 × 10,000 pixels, and the data type is 16-bit signed integer, consistent with the input images. The Vegetation Index (VI) values are stored as integers in the range −10,000 to 10,000 and are interpreted in physical units by applying a VI scale factor of 0.0001, corresponding to a real NDVI range of −1.0 to 1.0. The Cloud image uses binary values, where 1 indicates cloud and 0 indicates clear. The Flag image encodes, in MMDD format, the date on which the representative NDVI value was selected for each pixel. Within each image, missing pixels and ocean pixels are fixed to special constants so that they can be consistently identified in downstream processing: the Nodata value is set to −32,768, and the Sea value is set to−10,001.

2.4.3. Metadata Specification

A separate Meta file stores the quality indicators and ancillary information produced during the compositing and gap-filling process. Major quality indicators include the number of valid input images used during the composite period (“COMPOSITE_IMAGE_COUNT”), the residual cloud ratio after periodic compositing (“RESIDUAL_CLOUDPIXEL”), the ratio of valid pixels within the tile, and the ratio of pixels for which the period-representative image was used.

In addition, the metadata file records general product information such as product level (L3VP), tile ID, coordinate reference system, and spatial resolution, composite period, algorithm version used, and the definitions of the Nodata and Sea values. This information allows the operational system and subsequent application modules to automatically interpret and use the composite products (Table 1).

Table 1. Metadata table for periodic composite vegetation index products

2.4.4. Implementation and Runtime Performance

From an implementation perspective, the periodic compositing and gap-filling modules were developed in a Python 3.x environment. The system uses NumPy and SciPy for numerical computation, scikit-learn for training and applying the GPR model, and GDAL/rasterio for large-volume GeoTIFF input/output and coordinate transformation.

The periodic compositing module is designed to interface with the operational system of the agricultural satellite through an Application Order (AO) mechanism. When an AO file in XML format specifies the input paths (L3V NDVI and cloud), output paths, composite period, tile IDs, and other necessary parameters, the compositing process is executed automatically in batch mode.

Based on performance measurements obtained during trial operations in 2023, the time required to generate a 10 day composite set for a single tile on a typical server environment is approximately 5 minutes, which is sufficient for near-real-time operation considering the satellite’s revisit cycle and processing schedule.

The spatial and temporal gap-filling modules are not mounted directly in the core operational chain but are implemented as independent stand-alone modules that can be executed by users who wish to fill gaps in the periodic composite NDVI. The time required to perform the spatial and temporal corrections for one tile is also on the order of 5 minutes. Owing to these characteristics, the proposed framework can be used as a flexible NDVI production system from the early test-operation phase of the agricultural satellite and can be adapted to future changes in sensor specifications or observation strategies.

3. Results

3.1. Visual Results of Periodic Compositing and Spatio-Temporal Gap-Filling

We applied the proposed periodic compositing and spatio-temporal gap-filling algorithm to five standard tiles (Tile IDs: 96, 97, 109, 121, 130) from July to September 2023 to evaluate its performance. For each tile, we used three daily NDVI images with different cloud conditions within a 15 day compositing window as input. This 15 day interval is specific to the tile-based case study shown in this section and is used to illustrate the behavior of the CC-MNC and gap-filling modules under different cloud conditions; the operational NDVI time series adopts a fixed 10 day compositing window. The acquisition dates of the selected NDVI scenes are summarised in Table 2: 2023-07-06, 2023-07-21, and 2023-08-05 for Tiles 96, 97, and 109; and 2023-08-25, 2023-09-09, and 2023-09-24 for Tiles 121 and 130. The cloud fraction of the individual scenes ranged from 0.0% to 27.96%. Among the five tiles, Tile 130 had the highest combined cloud fraction (32.40%), so it was chosen to visually examine the performance of the compositing and gap-filling procedures (Table 3).

Table 2. Cloud fraction of three daily NDVI scenes for each tile

Table 3. Cloud amount for three daily NDVI images of tile 130

Fig. 7 shows the three daily NDVI images for Tile 130 on 25 August 2023 (Fig. 7a), 9 September 2023 (Fig. 7b), and 24 September 2023 (Fig. 7c). Despite corresponding to the late growing season after the monsoon, all three scenes contain extensive cloud cover, resulting in discontinuous NDVI patterns over paddy and upland fields. The corresponding cloud masks (Fig. 8) and the union of the cloud-covered areas are summarised in Table 3. The cloud fractions of the three scenes are 24.88%, 3.74%, and 5.60%, respectively, and the union of the three masks yields 28,668,004 cloudy pixels (32.40%) out of 88,494,364 pixels over Tile 130.

Fig. 7. The daily NDVI of tile 130: (a) 2023.08.25, (b) 2023.09.09, and (c) 2023.09.24.

Fig. 8. Cloud area of tile 130: (a) 2023.08.25, (b) 2023.09.09, and (c) 2023.09.24. Blue shaded areas indicate cloud pixels.

This indicates that time-series analysis of agricultural areas based on single-date NDVI images alone would be severely constrained, even with a 15 day compositing interval. Fig. 9 shows the same tile and presents the following: (a) the CC-MNC-based periodic composite NDVI and residual clouds, (b) the result after the first stage of spatial gap-filling and the remaining clouds, (c) the final NDVI and the residual clouds after the second stage of temporal gap-filling. The composite result (Fig. 9a) shows that cloud-free observations from the three dates are selected such that large-scale cloud bands are substantially removed, and field-scale NDVI patterns become more spatially continuous.

Fig. 9. NDVI compositing and spatio-temporal gap-filling results with residual clouds for tile 130.

However, small patch-like gaps still remain in areas where clouds persist throughout the entire compositing period or where the viewing-angle constraint is not satisfied. The first-stage spatial gap-filling result (Fig. 9b) fills gaps of ≤ 5 pixels using the median NDVI of surrounding cloud-free pixels, which is more robust to outliers such as residual cloud-contaminated pixels. This smooths irregular, speckled gap patterns along cloud boundaries while preserving the overall field structure. Finally, the second-stage temporal gap-filling result (Fig. 9c) corrects medium- to large-scale missing regions by exploiting NDVI time series from preceding and subsequent composite periods. As a result, an NDVI field that is continuous across the entire tile and free from cloud-induced gaps is obtained. Throughout this process, fine-scale spatial structures, such as paddy–upland boundaries, rivers, irrigation channels and forest–cropland interfaces, are preserved while only cloud-affected areas are reconstructed naturally.

3.2. Definition of Gap-Recovery Rate and Evaluation Methodology

To quantitatively evaluate the performance of the proposed periodic compositing and spatio-temporal gap-filling algorithm, we used the three daily NDVI images per tile described in Section 3.1 for the five standard tiles (Tile IDs: 96, 97, 109, 121, 130) shown in Fig. 1. The cloud fractions of NDVI images #1, #2, and #3 range from 0.0% to 27.96%, while the cloud fraction of the integrated missing area—defined as the union of the three cloud masks—ranges from 13.61% to 32.40% across tiles, with an average of 21.37% (Table 2). This integrated cloud mask is treated as the “initial” missing area for the compositing period.

In this study, we define two complementary indicators to assess the algorithm: (1) the compositing-stage gap-recovery rate, which measures how effectively the CC-MNC compositing step alone reduces cloud-induced gaps; and (2) the correction-stage gap-recovery rate, which measures how effectively the subsequent spatial and temporal gap-filling procedures recover the residual gaps left after compositing.

First, the compositing-stage gap-recovery rate is defined as follows. Let N_init denote the number of cloudy pixels in the integrated cloud mask (i.e., the total initial missing area within the compositing period) and N_{res, comp} denote the number of cloudy pixels remaining in the residual cloud mask after applying the CC-MNC periodic compositing. Then, the compositing-stage gap-recovery rate R_comp is given by

\(\begin{align}R_{\text {comp }}(\%)=\left(1-\frac{N_{\text {res, comp }}}{N_{\text {init }}}\right) \times 100\end{align}\)       (2)

Here, R_comp represents the proportion of the initial cloud-induced missing area that is recovered during the periodic compositing step; values closer to 100% indicate that the compositing process itself has effectively mitigated cloud-related gaps. Next, the correction-stage gap-recovery rate is defined by treating the residual cloud mask after CC-MNC compositing as the new “total” missing area to be corrected by the spatial and temporal gap-filling modules. Specifically, let N_{res, comp} (as above) denote the number of cloudy pixels remaining after the compositing step, and N_unfilled denote the number of pixels that still remain cloudy (i.e., uncorrected) even after applying both the spatial gap-filling and the temporal GPR-based gap-filling procedures.

Then, the correction-stage gap-recovery rate R_corr is defined as

\(\begin{align}R_{\text {corr }}(\%)=\left(1-\frac{N_{\text {unfilled }}}{N_{\text {res, comp }}}\right) \times 100\end{align}\)       (3)

In this formulation, R_corr quantifies how effectively the combined spatial and temporal gap-filling stages recover the residual missing area left after periodic compositing. A value of 100% indicates that all clouds remaining after CC-MNC compositing have been successfully removed by the correction stages, whereas lower values imply that some residual gaps remain. Together, R_comp and R_corr provide a clear stage-wise assessment of the effectiveness of the compositing and subsequent correction processes.

3.3. Quantitative Performance of Periodic Compositing and Gap-Filling

Based on the integrated cloud masks defined in Section 3.2, the initial missing area for each tile was derived from the union of the three daily NDVI cloud masks from July to September 2023. As summarised in Table 4, the initial cloud cover ranged from 13.61% (Tile 109) to 32.40% (Tile 130), averaging 21.37% across the five tiles. This indicates that more than one-fifth of the agricultural area within the tiles was obscured by clouds on average during the evaluation period, highlighting the need for effective gap reduction. After applying CC-MNC-based periodic compositing, the average residual cloud cover decreased to 0.93% (Table 4). Converting this to the gap-recovery rate relative to the initial integrated missing area shows that the compositing step achieved an average recovery rate of 95.50%.

Table 4. Cloud cover and gap-recovery rates before and after compositing and gap-filling for each tile

This demonstrates that compositing alone can remove most cloud-induced gaps. At the tile level, cloud cover in Tile 109 decreased from 13.61% to 0.04% (corresponding to a recovery rate of 99.71%), and Tile 130, which had the largest initial cloud cover, at 32.40%, achieved a recovery rate of 94.50%, reducing residual cloud cover to 1.78%. Tiles 96 and 97 also showed a significant improvement, with cloud cover decreasing from 15.62% and 31.61% to 0.21% and 0.96%, respectively. Applying the same spatial gap-filling (first-stage correction) reduced the average residual cloud cover across the five tiles further to 0.47% (Table 4). The average gap recovery rate of this first correction step, computed relative to the residual gaps after compositing, was 97.54%, indicating that most of the remaining gaps were successfully filled using spatial information.

However, small residual gaps persisted in areas of complex terrain and land cover transitions, such as coastal zones and forest–cropland boundaries. After the second-stage temporal gap-filling using GPR, the residual cloud cover in all tiles reached 0%, yielding a final gap-recovery rate of 100% for every tile (Table 5). In other words, a fully gap-free NDVI time series was obtained for the agricultural areas once both correction stages had been applied, with no remaining cloud-induced missing pixels. Even in the unfavourable case of Tile 130, where initial cloud cover exceeded 30%, the residual cloud fraction was reduced to 0% in the final product. This confirms the high practical applicability of the proposed framework for operational use from qualitative and quantitative perspectives.

4. Conclusions

In this study, we developed an integrated multi-stage framework to produce reliable, cloud-free NDVI time series for agriculture by combining a CC-MNC periodic compositing approach with subsequent spatial and temporal gap-filling techniques. The proposed algorithm was tested on five representative agricultural grid tiles in South Korea (Tile IDs: 96, 97, 109, 121, 130) during the cloudy 2023 monsoon season. The results demonstrate a dramatic reduction in cloud-induced gaps at each stage of the process. Starting from an initial average cloud cover of ~21.37% in the input Sentinel-2 scenes, the CC-MNC compositing alone reduced the average missing-data fraction to only 0.93%.

This corresponds to about 95.5% of the cloud-induced gaps being filled in the first stage. After applying the spatial gap-filling step, the average remaining gap was further lowered to 0.47%, achieving a 97.5% total gap recovery up to that point. Finally, the GPR-based temporal gap-filling stage eliminated all remaining gaps, resulting in 0% residual cloud cover in the NDVI composites. Even the cloudiest case, tile 130 with an initial 32.40% cloud cover, was completely gap-filled to 0% missing data, confirming both visually and quantitatively that the framework can effectively resolve severe cloud-related data loss.

In addition to its effectiveness, the framework is computationally efficient. Generating a 10 day NDVI composite for one tile takes on the order of only 10 minutes on typical server hardware. This speed is well within the 3 day revisit interval of the upcoming CAS500-4 agriculture satellite, indicating the approach can meet near-real-time operational requirements. In fact, by using Sentinel-2 data as a proxy for CAS500-4, we verified that the system can be deployed on the CAS500-4 standard grid from the very start of the satellite’s operations.

From a scientific perspective, this study demonstrates that integrating CC-MNC compositing with local spatial and GPR-based temporal gap-filling can significantly improve gap recovery and temporal continuity compared with conventional MaxNDVI-only or purely temporal interpolation approaches, while preserving field-level spatial patterns in heterogeneous Korean croplands. Operationally, the proposed framework provides a practical blueprint for implementing CAS500-4 NDVI production on the C4SW Standard Grid in near-real time, thereby supporting national-scale crop monitoring and decision-making under persistently cloudy conditions.

In summary, the proposed NDVI time-series production system can reliably provide complete gap-free NDVI maps even during cloud-prone growing seasons. It is expected to serve as a practical remote-sensing tool for crop monitoring in operational programs, including the new CAS500-4 satellite mission, thereby contributing to improved agricultural monitoring under persistently cloudy conditions.

Author Contributions

Conceptualization: Kim JY, Kim SH, Kim TH; Data curation: Kim JY, Kim SH, Eun J; Formal analysis: Kim SH, Eun J, Yoon SJ, Ryu HS, Kim TH; Funding acquisition: Kim TH; Investigation: Yoon SJ, Ryu HS, Kim TH; Methodology: Kim SH, Eun J, Kim TH; Project administration, Supervision: Kim TH; Validation: All authors; Writing–original draft: Kim JY, Kim TH; Writing– review & editing: All authors.

Conflicts of Interest

No potential conflict of interest relevant to this article was reported.

Funding

This work was carried out with the support of the “Cooperative Research Program for Agriculture Science and Technology Development (Project No. RS-2021-RD009991)” by the Rural Development Administration, Republic of Korea.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.