Landsat and Sentinel-2 data archives provide ever-increasing amounts of satellite data. However, the availability of usable observations greatly varies spatially and temporally. Pixel-based compositing that generates temporally equidistant cloud-free synthetic images can mitigate temporal variability, by constructing uninterrupted time series using different compositing windows. Here, we evaluated the feasibility of using compositing windows ranging from five days to one year for 1984-2021 Landsat and 2015-2021 Sentinel 2 time series to derive uninterrupted time series across Europe. We considered separate and joint use of both data archives and analyzed the spatio-temporal availability of composites during each calendar year and pixel-specific growing season across a variety of time windows and hypothesizing data interpolation. Our results demonstrated opportunities and limitations in the available data records to support medium- and long-term analyses requiring uninterrupted time series of composites with sub-annual temporal resolution. Spatial disparities across different compositing windows provide guidance on the feasibility of workflows relying on different data densities and on the challenges in wall-to-wall analyses. The feasibility of consistent time series based on composites with sub-monthly aggregation periods was mostly limited to the combined Landsat and Sentinel-2 archives after 2015, yet in some geographies requires interpolation of up to 50% of data.

We used all Landsat surface reflectance Level 2, Tier 1 (Collection 2) scenes from 1984 through 2021 and Sentinel-2 TOA reflectance Level-1C (pre‑Collection-1; European Space Agency, 2021) scenes from 2016 through 2021 acquired over Europe, as available in Google Earth Engine (data accessed in June 2022; Gorelick et al., 2017). We utilized Seninel-2 Level-1C data instead of Level-2A because the Level-2A inherent quality data lack the desired scope and accuracy (Baetens et al., 2019; Coluzzi et al., 2018). Yet, the Level-1C products are accompanied by cloud probabilities (Zupanc, 2017) facilitating improved cloud screening. Furthermore, for cloud screening we also used Band 10 (Cirrus), which is not available as Level‑2A. Because we performed data availability analyses, i.e., we tallied the daily presence/absence of usable observation, the disparity between Landsat and Sentinel-2 reflectance values was here irrelevant, and the intra-sensor normalization was not needed. The difference in processing levels, however, played out in cloud, shadow, and snow masking accuracy, where the Sentinel-2 workflow assembles several approaches with known accuracies (Skakun et al., 2022), but has not been evaluated as a whole. We acknowledge that for real-life reflectance-based applications, data from corresponding processing levels need to be used and the reflectance normalized among the sensors (Okujeni et al., 2024). We recommend thus either preprocessing of the Sentinel-2 TOA data to achieve the desired quality of masks, or linking Sentinel-2 Level-A2 data with Level-1C band 10 and relevant Cloud Probability scenes for more rigorous cloud screening.

To ensure that only pixels with the highest quality entered the analysis we applied conservative pixel-quality screening. For Landsat scenes, we excluded all pixels flagged as cloud, shadow, or snow using the inherent pixel quality bands (Foga et al., 2017; Z. Zhu & Woodcock, 2012) and discarded saturated pixels (Zhang et al., 2022). We further used the quality bands to exclude all data gaps in the Landsat 7 acquisitions occurring due to the SLC scanline failure (Andréfouët et al., 2003). Although the accuracy of the inherent pixel‑quality bands differs among the Landsat sensors due to the differences in the sensor’s build and thus availability of thermal and cirrus‑specific bands (Foga et al., 2017), the Landsat quality bands are acclaimed standardized quality product. Finally, owing to Landsat 7’s orbit drift (Qiu et al., 2021), we excluded all ETM+ scenes acquired after 31^st December 2020.

We used a 20‑km grid of 16,642 equidistant points to analyze the availability of useable Landsat and Sentinel-2 observations over Europe. We distributed points according to the Lambert azimuthal equal-area projection (LAEA, EPSG:3035), which is the preferred projection for EU-wide products. Despite LAEA being the equal-area projection, the distance distortion within our study area was mostly below 10 m, which is less than one pixel in high‑resolution Sentinel-2 bands. The systematic gridded sampling design ensured good representation of the West-East and South‑North climatic and phenological gradients, and facilitated graphical presentation of results.

We derived the time series of usable Landsat and Sentinel-2 observations over Europe sampling individual pixels spaced systematically every 20 km in the latitudinal and longitudinal directions. We identified sampling locations according to the Lambert azimuthal equal-area projection (LAEA, EPSG:3035), which is the preferred projection for EU‑wide products. Despite LAEA being the equal-area projection, the distance distortion within our study area was mostly below 10 m, which is less than one pixel in high‑resolution Sentinel-2 bands. The systematic point sampling design is used to derive overview statistics for big datasets and in nearest neighbor-based rescaling of rasters. The 20-km sampling interval resulted in 16,642 locations over land ensuring good representation of the West-East and South‑North climatic and phenological gradients, as well as facilitating graphical presentation of results.

For each sampled pixel we recorded the date of the valid cloud-, shadow, and snow-free Landsat and Sentinel-2 acquisition. We used the information at the original resolution and assumed each sampled pixel to be a probabilistic sample of the surrounding 20x20-km area, making the process analogous to the nearest neighbor resampling. We excluded duplicated data entries coming from the vertical overlaps among Landsat tiles in the same row, and vertical and horizontal overlaps among Sentinel-2 granules from the same swath. This resulted in daily data availability for 1984-2021 (1 – valid observation; 0 – no data or no valid observation), which we used to derive availability information for composites with aggregation periods of five, 10, 15, 20, and 25 days; one, two, three, four, six and 12 months. The non-overlapping compositing windows compartmentalized daily information for each year into 73, 37, 24, 18, 15, 12, six, four, three, two, and one composites for each calendar year, respectively. We used January 1^st as the starting date for the compositing window sequence for each year. When the last compositing window was shorter than half its window width, we merged it with the penultimate composite. For each data point and every considered aggregation period we recorded the amount of available observations and considered a composite as ‘successful’ if at least one valid observation was available.

The data are distributed as a csv and GeoTIFF (both formats comprising exactly the same information) and can be open and query using any software able to handle these data formats.