Continuous, consistent, and long time-series from remote sensing are essential to monitoring changes on Earth's surface. However, analyzing such data sets is often challenging due to missing values introduced by cloud cover, missing orbits, sensor geometry artifacts, and so on. We propose a new and accurate spatio-temporal prediction method to replace missing values in remote sensing data sets. The method exploits the spatial coherence and temporal seasonal regularity that are inherent in many data sets. The key parts of the method are: 1) the adaptively chosen spatio-temporal subsets around missing values; 2) the ranking of images within the subsets based on a scoring algorithm; 3) the estimation of empirical quantiles characterizing the missing values; and 4) the prediction of missing values through quantile regression. One advantage of quantile regression is the robustness to outliers, which enables more accurate parameter retrieval in the analysis of remote sensing data sets. In addition, we provide bootstrap-based quantification of prediction uncertainties. The proposed prediction method was applied to a Normalized Difference Vegetation Index data set from the Moderate Resolution Imaging Spectroradiometer and assessed with realistic test data sets featuring between 20% and 50% missing values. Validation against established methods showed that the proposed method has a good performance in terms of the root-mean-squared prediction error and significantly outperforms its competitors. This paper is accompanied by the open-source R package gapfill, which provides a flexible, fast, and ready-to-use implementation of the method.