Data selection
Selecting which observations to include requires navigating two interrelated concerns: the temporal distribution of acquisitions (Temporal distribution) and the quality criteria for both optical (Optical data selection) and SAR data (SAR data selection). They involve trade-offs that are straightforward to state but difficult to resolve in practice, particularly at global scale.
Temporal distribution
We frame the present discussion around the task of producing a yearly map. Under this setting, we discuss how to temporally select EO data within a period of approximately 12 months. There is no one-size-fits-all guideline for temporal distribution, much of it ultimately depends on how much data one can afford to process. When deciding how many time steps to obtain, three broad strategies emerge: single time-steps, time-series, or composites, each with its own trade-offs. Time-series are arguably the richest option, capturing phenomena like seasonal variation across a defined period. [1] showed that a monthly Sentinel-2 time-series of one-year preserved the seasonal nuances and geolocation shifts that median composites tend to wash away, and [2] found benefits in using a time-series of composited PlanetScope imagery over using a single image. The trade-off is data size: they are memory-intensive to store by nature. Furthermore, for each timestep in the time series comes the question of which observation to choose. In the case of optical data, most researchers opt for a time-series with a monthly temporal resolution, choosing the least cloudy product per month [1], [3], [4], [5]. Composites offer a pragmatic alternative, distilling an entire period of observations into the memory footprint of a single time-step, as in [6], [7], [8], [9], [10], [11]. Yet, compositing is not without its own costs: it requires deliberate design choices and non-trivial computations. Typical approaches consist of performing mean or median compositing, the latter involving computing percentiles, and the former being more sensitive to outliers. Moreover, if compositing cannot be performed on-the-fly before downloading the data, practitioners must store (at least temporarily) the time-series the composite will be computed from, undermining much of the storage advantage. While the per-pixel processing cost may appear modest in isolation, it scales dramatically: at a global scale, every additional step must be applied to each pixel on Earth, which can quickly become prohibitive. One can also produce a time series of composites, e.g., monthly composites of Sentinel-2 [4], seasonal composites [12], [13], [14], [15], or using a predetermined moving window [2], [16], [17], [18]. Single time-steps, by contrast, are the most cost-efficient route but they come at the expense of the temporal dimension, which often carries rich, discriminative information. Indeed, [18] found that performance steadily improved with access to longer input timespans when mapping tree mortality, and [19] found similar results for canopy height mapping.
Beyond the choice of temporal representation, the time range itself must be defined. We focus on vegetation-related tasks throughout this section, where the trade-offs are particularly well documented, but the same considerations apply across domains. A common default is to restrict the analysis to the leaf-on season [5], [6], [18], [20], [21], [22], but identifying that period for every pixel on Earth is in itself non-trivial. Operational phenology products offer one solution: [18] use the VIIRS Global Land Surface Phenology product [23], while [21] use the MODIS equivalent [24]. Both products are released at annual intervals, so inter-annual shifts in phenological dates are in principle captured. However, the underlying algorithms assume a vegetation cycle with at most two growing cycles per year, an assumption that breaks down for multi-cropping systems with more than two cycles. [15] bypass this limitation by relying on global crop calendars instead [25]. More fundamentally, the annual leaf-on/leaf-off framing is itself a temperate-zone assumption that does not generalize to tropical systems [26]. The leaf-on framing can also be limiting in its own right: [1] show that incorporating both summer leaf-on and winter leaf-off observations improves performance, and [27] find that for biomass estimation in Finland, early-spring Sentinel-2 acquisitions are more informative than summer ones, which they attribute to “the higher contrast between vegetation and snow, which enhances the ability to distinguish vegetation properties”. This logic does not extend everywhere, but regions with strong seasonal contrast can benefit from year-round data. A more task-specific variant samples contrasting images from different times within the same year to expose intra-annual variation, as applied by [28] for crop field segmentation. The calendar year itself is also a poor natural boundary: in the Southern Hemisphere, summer spans the end of one year and the beginning of the next, so a global-scale pipeline must accommodate such cross-year shifts rather than assume that a single calendar year is a sufficient window.
In summary, the temporal strategy should be guided by the specifics of the task at hand, keeping in mind that any processing will have to be applied to every pixel on Earth and that data volume scales with the size of the land mass. Beyond when and how often to observe, the next question is which individual observations are reliable enough to include, a decision governed by the quality criteria discussed next.
Optical data selection
Even after temporal selection, individual observations vary in usability. Two factors dominate optical imagery: cloud contamination and sensor-specific artifacts that introduce gaps or duplications. Rather than discussing the per-pixel quality flags shipped with each product (e.g., the Sentinel-2 SCL band or Landsat QA band), we focus on higher-level quality criteria that govern whether entire products should be included or discarded.
Upon selecting which data product to work with, a common filter in optical imagery is the cloud cover. We elaborate on the various methods to perform cloud detection in Section Data Preprocessing and henceforth assume access to a reliable cloud mask. A common approach is to discard data products with a cloud cover above a certain threshold (e.g., 10% in [6], 35% in [29], 60% in [12], 75% in [30], 90% in [28]). However, for regions with persistent cloud cover, filtering out entire data products based on their cloud cover percentage induces the risk of throwing away valuable non-cloudy pixels. For example, when compositing, products with high cloud cover may contain cloud-free areas that fill critical gaps in the composite. There is no straightforward way to best select pixels based on cloud masks, and practitioners have developed their own heuristics based on the task at hand. For example, [30] designed a Sentinel-2 product selection algorithm that maximizes spatial coverage while minimizing redundancy when creating cloud-free composites: products are chosen based on their contribution of valid pixels to under-represented spatial regions within a tile. Computing or accessing a cloud mask at the stage of data product selection can be computationally expensive, since it must be done for each candidate product. Instead, most optical data providers report cloud cover statistics in their products metadata. For most sensors, this cloud cover percentage is computed over the actual data footprint, either because their products follow the orbital swath (Landsat 8/9, Sentinel-3, Planet, Airbus) or because they explicitly exclude undetermined pixels from the calculation (MODIS, VIIRS). However, particular care must be taken when working with Sentinel-2 data, where the CLOUDY_PIXEL_PERCENTAGE is computed over the entire 110km \(\times\) 110km tile, which may include NODATA regions and lead to artificially low cloud cover values for partially covered tiles.
Another important criterion for selecting optical data products is the occurrence of artifacts. For example, Landsat 7 imagery acquired after 2003 exhibits wedge-shaped gaps from the Scan Line Corrector failure, affecting roughly 22% of each scene [31], [32]. In what follows, we focus on Sentinel-2 data products, as they are among the most widely used sources of freely available optical data. Sentinel-2 artifacts manifest mainly in two ways: incomplete orbits and split partial products. Because the distribution unit (granule) is not geometrically aligned with the sampling unit (orbit), each tile is covered by multiple orbits, some of which encompass it entirely. This is illustrated in Figure 4.1 (top), where tile 31TGN is partially covered by orbit 008 and fully covered by orbit 108. In such cases, the appropriate orbit can be identified directly from the MTD_MSIL2A.xml file of each product: the NODATA_PIXEL_PERCENTAGE should be 0.0 for the orbit that fully covers the tile, and >0.0 for the other. Most tiles, however, are not fully covered by any single orbit. These tiles require mosaicking data from multiple orbits, which may have been acquired at different time steps. Tile 30SUH provides a typical example: as shown in Figure 4.1 (middle), two products acquired over two days are needed to obtain a complete image.
A more unusual case occurs when a tile that should be fully covered by a single orbit yields two products sharing the same date, orbit, and footprint. One cause is a change in processing baseline, in which case the product with the highest baseline can simply be retained. When the baselines are identical, however, the duplication indicates a data strip split: a single continuous observation broken into separate segments during ground processing, slicing a tile into two partial files that must be mosaicked to reconstruct the full area (Figure 4.1, bottom).
Addressing these issues is important to avoid duplication of pixel values, which would bias any statistics computed over the data — whether during compositing or when calculating overall dataset statistics from, for example, all Sentinel-2 products acquired within a given time frame, as in [33], [34]. To mitigate this, [16] capped acquisition to at most one observation per day per tile, retaining the observation with the highest number of valid pixels when multiple were available. While sensible, this strategy addresses partial orbits only in part: it does not account for split partial products, nor does it guarantee de-duplication of pixels, an issue closely tied to the broader question of EO data redundancy [35], which we discuss in Section Spatial gridding.
SAR data selection
SAR imagery is not affected by clouds, but its side-looking acquisition geometry introduces its own artifacts: geometric distortion due to foreshortening, layover, and shadow (see Figure 4.2). These deterministic perturbations arise from the interaction between terrain relief and the radar’s look vector: given a digital elevation model and the acquisition geometry, the distortion pattern can be predicted.
For Sentinel-1 in Interferometric Wide (IW) mode, three factors define this viewing geometry: the relative orbit number or track (the satellite’s repeating ground path), the sub-swath (which sets the incidence angle range), and the pass direction (ascending or descending). Ascending and descending acquisitions of the same area are complementary (slopes that suffer layover under one geometry are often well-imaged under the other) and combining the two is a standard strategy for filling distortion gaps in mountainous regions. More generally, several tracks are typically available over any target area, each contributing hundreds of acquisitions over the mission’s lifetime. Both the data volume and the need for geometric consistency over time (required for analyses such as interferometric stacks, amplitude change detection, or feature extraction) may push toward the selection of a single track. ESA’s Sentinel-1 Burst ID Map [36] provides the metadata needed to make this choice. It consists of global polygons that divide each sub-swath into smaller units called bursts, each annotated with its track and sub-swath identifier. Together with the pass direction, these identifiers specify the viewing geometry for any given location, and three approaches of increasing detail can be used to select among the candidates. The simplest uses the burst grid alone: candidates are filtered by sub-swath, typically with a preference for mid-range incidence angles, which is sufficient for flat or moderate terrain. For rougher terrain, combining the burst grid with a digital elevation model allows candidates to be ranked by the alignment between the radar’s look direction and the dominant slope aspect. For a more exact comparison, tools such as SNAP’s SAR Simulation Operator can simulate the per-pixel layover and shadow masks that would arise from each candidate, though the per-candidate simulation cost makes this approach impractical at scale.
Once the temporal strategy and quality criteria for data selection have been established, the next challenge is to process the selected observations into a consistent, analysis-ready format.