Predicting missing values in spatio-temporal satellite data

1. Introduction

• Atmospheric disturbances and off-nadir viewing can introduce measurement errors or omissions in remotely sensed data.Remote sensing is used to study Earth-surface processes and offers broad spatial and temporal coverage compared with ground-based measurements.Source: §1 · p. 3

1.1. Missing values in satellite data

• Good-quality MODIS data can have substantial gaps, reaching 100% of pixels in certain Alaska-region images.The introduction notes that such missing data may negatively influence inference about Earth-surface processes or render some analyses impracticable.Source: §1.1 · p. 3

1.2. Existing approaches to handle missing values

• Missing-data strategies either use analyses robust to gaps or predict missing values before analysis.Prediction can be applied to a data product with missing values or integrated into preprocessing.Source: §1.1 · p. 4, §1.2 · p. 4

1.2.1. Robust analysis methods

• Vegetation time-series methods include smoothing, Fourier and LOESS approaches, and specialized breakpoint algorithms; comparative studies report different best-performing methods across scenarios.Bolin et al. (2009) use both temporal and spatial correlation to study spatial patterns of temporal trends in NDVI data.Source: §1.2.1 · pp. 4–5

1.2.2. Reconstructing a complete data product

• Reconstruction methods range from spatial Kriging and Co-Kriging to spatio-temporal models and non-geostatistical pixel-based or regression approaches.Some regression approaches use additional land-cover data or MODIS quality information.Source: §1.2.2 · pp. 5–6

1.3. Outline

• The paper formalizes a generic spatio-temporal subset-prediction procedure, proposes an accurate instance with uncertainty quantification, and outlines evaluation across four scenarios against TIMESAT and Gapfill-Python.The proposed procedure reconstructs a complete data product before satellite-data analysis.Source: §1.3 · p. 6

2.2. Formal description

• Each missing value is predicted separately from a four-dimensional array indexed by spatial location, seasonal position, and year.The procedure can be applied to a wide range of remotely sensed data.Source: §2.1 · p. 7

2.1. Illustration of the concept

• The subset-predict procedure handles values one at a time, using only a local subset to support scalable, distributable computation.Its adaptive subset can recover local features such as shifted seasonal patterns or long-term trends.Source: §2.1 · p. 9

2.2.1. Subset to a neighborhood

• The neighborhood expands spatially until it contains enough non-empty images and enough observed values in the target image.The initial neighborhood also spans specified seasonal and yearly ranges; the tuning parameters can be adjusted for a dataset.Source: §2.1 · p. 9, §2.2.2 · p. 12

2.2.2. Predict based on a subset

• Linear quantile regression predicts the missing value from an intercept and image rank, with the rank serving as a non-parametric seasonality adjustment.The quantile is estimated from images observed at the target location, using spatial neighbors if too few such values are available.Source: §2.1 · p. 10, §2.2.2 · p. 14

2.1. Illustration of the concept

• The open-source R package gapfill has a flexible design for tailoring the procedure and running it in parallel through openMP and MPI back-ends.Users can also construct new subset-predict algorithms with little effort.Source: §2.1 · p. 10

2.2.2. Predict based on a subset

• Subset images are ranked through pixel-wise comparisons of their non-missing values, assuming similar but potentially shifted spatial distributions.The scoring is designed to be only marginally affected by missing values.Source: §2.1 · p. 10, §2.2.2 · p. 13

2.3. Prediction uncertainties

• Uncertainty is assessed for four procedure steps, but the combined prediction interval includes only ranking and quantile-estimation uncertainty because interactions make combining all four nontrivial.The combined interval reflects local spatial and temporal heterogeneity; resampling-based approaches can be inaccurate when their assumptions are unmet, as with missing-at-random for these data.Source: §2.3 · pp. 15–16

3.1. Data

• MODIS NDVI data from 2004–2009 were restricted to the seasonal period from day 145 to day 257, when each considered day had at least 30% good-data values.The study region is in northern Alaska, where winter has little available data because of missing sunlight and snow cover.Source: §3.1 · p. 17

3.2. Test scenarios

• Four 100 × 100-pixel validation scenarios retained the temporal structure and used realistic missingness patterns to reach 20%, 30%, 40%, and 50% missing values.Values were artificially removed from a region with about 12% missing data, leaving most true values known for comparison.Source: §3.2 · pp. 18–19

3.3.1. Predictions accuracy

• Prediction accuracy was evaluated using the proportion of values left missing, visual inspection, and RMSE against known values in the test scenarios.RMSE is calculated over non-missing predictions and compares predicted values with true values.Source: §3.3.1 · pp. 19–20

3.3.2. Uncertainty assessment

• Uncertainty components were compared by prediction-interval length, and the combined 90% interval was assessed for spatial and temporal patterns and coverage.Component analysis used days 145 and 161 from the 40%-missing scenario; coverage was evaluated on the full 40%-missing scenario.Source: §3.3.2 · p. 20

3.3.3. Comparison with TIMESAT and Gapfill-Python

• The comparison evaluated gapfill against Gapfill-Python and TIMESAT, with TIMESAT configured to use a double-logistic smoother.The authors note that results could improve with other tuning settings, so the comparisons provide an impression of performance.Source: §3.3.3 · pp. 20–21

4.1. Predictions accuracy

• All missing values were predicted across the test scenarios without visually introduced spatial patterns, and small-scale NDVI features and seasonal variation were recovered.For the 40%-missing scenario, Figure 3 shows predicted NDVI for day 177 of 2006, with observed pixels in black.Source: §4.1 · p. 22, Figure 3 · p. 22
• Early-season images had higher average RMSE, while spatial RMSE patterns resembled temporal data variability; early dates also tended to have more missing pixels.Figure 3’s middle panel shows average RMSE by date and its right panel shows the spatial distribution of pixel-wise average RMSE.Source: §4.1 · pp. 22–23, Figure 3 · p. 22
• Prediction errors widened as missingness increased; their medians stayed at zero with slight positive skew, and pixels below NDVI 0.5 had lower accuracy.Figure 4 shows predicted-versus-true scatterplots above histograms of prediction differences for the four test scenarios.Source: §4.1 · p. 23, Figure 4 · p. 24
• Across the full northern Alaska region, gapfill filled all missing values, and visual inspection found no artificially introduced patterns.Source: §4.1 · p. 24

4.2. Uncertainty assessment

• The sorting step contributed the largest prediction-interval widths, followed by quantile estimation.Figure 5 summarizes uncertainty contributions from the procedure’s four main steps.Source: §4.2 · p. 24, Figure 5 · p. 25
• The prediction interval had 93% overall coverage, indicating slight average overestimation; its width followed a seasonal U-shape with deviations at the analyzed period’s boundaries.Figure 5 depicts average interval widths and daily coverage; the study considers days 145 through 257.Source: §4.2 · p. 25

4.3. Comparison with TIMESAT and Gapfill-Python

• gapfill had the lowest RMSE in all four scenarios; Gapfill-Python was comparable, while TIMESAT RMSEs were about two times higher.RMSE was calculated only for non-missing predictions; gapfill RMSE was also evaluated on the subsets available for the other methods.Source: §4.3 · p. 26, Table 2 · p. 26

• The hybrid method combines algorithmic selection and image scoring with quantile regression and permutation tests to pair speed and scalability with uncertainty-quantification tools.The discussion contrasts this design with classical geostatistical models whose computational workload may exceed available resources and algorithmic methods for which uncertainty quantification is difficult.Source: §5 · p. 27
• The algorithm handles each missing value separately, enabling workload distribution across computing units for large datasets.The implementation uses R’s foreach framework to interface with parallel back ends, including OpenMP and MPI.Source: §5 · p. 27
• Testing used observed MODIS NDVI data and observed missingness patterns to approximate real use cases; software comparisons remain relative to the chosen scenarios.Source: §5 · p. 28
• The open-source R package gapfill includes source code, documentation, and data examples; users can tailor the algorithm through parameters or custom subset and prediction functions.Source: §5 · p. 28

6. Conclusion

• The method handled a higher percentage of missing values and provided the most accurate predictions in terms of RMSE compared with two alternative gap-fill methods.Validation scenarios used MODIS NDVI data with 20%–50% missing values.Source: §6 · p. 28
• The open-source R package gapfill can be tailored to specific needs and is suitable for regularly spaced data products.The software works for the presented data and its flexible design allows users to adapt it with little effort.Source: §6 · p. 28
• Statistical components of the algorithm enable uncertainty quantification for predicted values.Source: §6 · p. 28
• Gap-filling improves or enables analyses that are designed for complete remotely sensed datasets.Source: §6 · p. 28