2.1. What are Earth System Data Cubes (ESDCs)?

The concept of ESDCs was introduced along with the Earth System Data Lab (ESDL, Mahecha et al., Reference Mahecha, Gans, Brandt, Christiansen, Cornell, Fomferra, Kraemer, Peters, Bodesheim and Camps-Valls2020), an integrated data and analytical hub that aimed to unify multiple heterogeneous Earth system data streams into a standard data model with a unique Coordinate Reference System (CRS). ESDCs represent multidimensional data structures designed to facilitate streamlined access, analysis, and manipulation of Earth system data. ESDCs comprise labels as dimensions defining the cube’s axes, an array of grids with their associated coordinate values distributed along these dimensions, and univariate data associated with each grid cell. Furthermore, in this paper, we add a new component: a suite of attributes that characterise the data, the dimensions, and the complete ESDC entity.

The dimensions are a set of labels describing the axes of the ESDC. Generally, these dimensions comprise space (e.g. “x” and “y”), time, and variables. Nevertheless, further dimensions can be added (e.g. “pressure levels”, “model ensembles” or “time series components”). It is crucial to emphasise that while ESDCs conventionally incorporate spatial and temporal dimensions (e.g. latitude, longitude, and time), they are not confined to this paradigm (cf. Table 1 of Mahecha et al., Reference Mahecha, Gans, Brandt, Christiansen, Cornell, Fomferra, Kraemer, Peters, Bodesheim and Camps-Valls2020). ESDCs can exhibit different dimensions, and the number of dimensions is called the order of the ESDC. Thus, an increment in the ESDC’s complexity according to its dimensions is given by their order (e.g. a spatio-temporal grid of a univariate ESDC has an order of 3, while the order of a multivariate ESDC is 4).

The grouping of grids consists of discrete subsets derived from the domain of each dimension’s axis. The values of these subsets are referred to as coordinates, and, in the case of a regular grid, they determine the data’s resolution along that specific dimension. For instance, a grid determining a resolution of 0.5 degrees for the “latitude” dimension in a global ESDC may have a set of coordinates $ \mathrm{grid}\left(\mathrm{latitude}\right)=\left\{-89.75,-89.25,\dots,, \hskip2pt ,\mathrm{89.25,},\hskip2pt ,89.75\right\} $ . While coordinates are often associated with numerical values (e.g. latitudes and longitudes), they can encompass a wide range of values. For instance, timestamps in a “time” dimension with a set of coordinates $ \mathrm{grid}\left(\mathrm{time}\right)=\left\{``2022-01-01",``2022-01-02",\dots, ``2022-12-31"\right\} $ , or components derived from a time series decomposition approach in a “component” dimension with a set of coordinates $ \mathrm{grid}\left(\mathrm{component}\right)=\left\{``\mathrm{raw}",``\mathrm{trend}",``\mathrm{seasonal}",``\mathrm{residual}"\right\} $ . The grids within an ESDC exhibit the following characteristics: 1) In the case of spatial dimensions, they reference the same CRS, 2) the coordinates within a grid share identical units, and 3) they must consist of at least two coordinates; otherwise, the dimension (and consequently the grid) is omitted. It’s important to note that, given these properties, irregular grids are also possible, with the temporal dimension grid being a typical example in EO data due to the irregular revisit times of some satellite missions (e.g. Sentinel-2).

The array of data represents scalar values corresponding to each grid cell. Typically, the data spans from observed measurements to modelled values. Nevertheless, one can also encounter higher-order features (that is data derived from operations performed on the original values), such as outcomes from time series decomposition or AI-generated products. Furthermore, flag values, which delineate data status, can be incorporated. Cells without data are denoted as “NA” (that is not available).

The collection of attributes comprises a series of key-value objects that provide additional details about the data. These objects serve as metadata and can offer descriptions ranging from individual variables (including their associated dimensions) to the entire ESDC. The information contained within these attributes typically encompasses a wide range of elements, including, but not limited to, names, acronyms, units, flag definitions, versions, and source details.

2.2. Relation of ESDCs to Image Collections and Data Cubes

Earth system data often exhibits heterogeneity and irregularity, particularly within EO data. Variability can manifest in different spatial resolutions, time units, projections, formats, and more, sometimes even within the same data product. Consequently, two robust approaches to retaining data integrity without succumbing to information loss due to transformations (e.g. reprojection, reduction, and resampling) are to utilise image collections (refer to Appel and Pebesma, Reference Appel and Pebesma2019 for a comprehensive distinction between image collections and conventional data cubes) or to adopt a process where original files are stored and indexed within a information-preserving data cube infrastructure based on their file metadata (Figure 1). In the latter approach, the original files can be stored locally or in the cloud while preserving the essential information intact.

Figure 1.

Representations of different storage systems for gridded data in Earth system research: Image collections (left), information-preserving data cubes (centre), and Earth system data cubes (ESDCs, right). Differences in these abstract representations have deep implications for data storage systems, accessibility, interoperability and metadata definitions.

A successful example of image collections is the GEE Catalogue. This extensive, multi-petabyte catalogue stores data in tiled images, where each image may encompass multiple bands, thereby preserving essential information. Furthermore, these images can be organised into an image collection if they share relevancy. GEE also offers the computational resources necessary for accessing and analysing their catalogued data. Within GEE, data cube-like operations can be seamlessly executed through dynamic on-the-fly reprojection, resampling, and reduction for the tiles where a subset of pixels was explicitly requested. It is worth noting, however, that users are required to conform to the specific Application Programming Interfaces (API) provided by GEE for processing and analysing the data effectively.

Standardising image collections and their access brings simplicity and promotes data usage across platforms. Currently, a widely recognised standard is the Spatio-Temporal Assets Catalog (STAC) specification. This specification empowers users to query data assets based on metadata and spatio-temporal criteria. Coupled with domain-specific API clients available for multiple programming languages (cf. Section 6.2) and GIS software (e.g. QGIS STAC PluginFootnote ⁶), users can easily retrieve data. The flexibility of the STAC specification has prompted numerous data providers to adopt it for creating their own data cataloguesFootnote ⁷, with notable examples including Microsoft Planetary Computer CatalogueFootnote ⁸ and the United States Geological Survey (USGS) Landsat Archive Catalogue (stored in the Amazon Simple Storage Service, S3)Footnote ⁹.

The information-preserving data cube approach is exemplified by the Open Data Cube (ODC) initiative, a prominent model in this field (Killough, Reference Killough2018; Killough et al., Reference Killough, Siqueira and Dyke2020)Footnote ¹⁰. This approach has played a pivotal role in informing governmental actions and policies, as evidenced by their integration into national and regional data cube frameworks (Dhu et al., Reference Dhu, Giuliani, Juárez, Kavvada, Killough, Merodio, Minchin and Ramage2019; Sudmanns et al., Reference Sudmanns, Augustin, Killough, Giuliani, Tiede, Leith, Yuan and Lewis2022). Noteworthy instances of these initiatives include Digital Earth Africa (DE Africa, formerly known as Africa Regional Data Cube, Killough, Reference Killough2019), Digital Earth Australia (DE Australia, previously Australian Geoscience Data Cube, Lewis et al., Reference Lewis, Oliver, Lymburner, Evans, Wyborn, Mueller, Raevksi, Hooke, Woodcock, Sixsmith, Wu, Tan, Li, Killough, Minchin, Roberts, Ayers, Bala, Dwyer, Dekker, Dhu, Hicks, Ip, Purss, Richards, Sagar, Trenham, Wang and Wang2017; Dhu et al., Reference Dhu, Dunn, Lewis, Lymburner, Mueller, Telfer, Lewis, McIntyre, Minchin and Phillips2017), the Colombian Data Cube (CDCol, Ariza-Porras et al., Reference Ariza-Porras, Bravo, Villamizar, Moreno, Castro, Galindo, Cabera, Valbuena and Lozano2017; Bravo et al., Reference Bravo, Castro, Moreno, Ariza-Porras, Galindo, Cabrera, Valbuena and Lozano-Rivera2017; Villamizar et al., Reference Villamizar, Castro, Ariza-Porras, Mancipe, Cabrera, Pachon, Ramirez, Fonseca, Lozano-Rivera, Cabrera and Becerra2018), and the Swiss Data Cube (SDC, Giuliani et al., Reference Giuliani, Chatenoux, Bono, Rodila, Richard, Allenbach, Dao and Peduzzi2017).

While both of these approaches excel in preserving data integrity and offering flexibility for various analyses, achieving Earth system data interoperability necessitates their integration into a unified structure through ESDCs. These ESDCs can be constructed from either approach. For instance, in the case of image collections, it is feasible to request pixels from GEE (Clinton, Reference Clinton2023), and data transformations can be executed within the GEE cloud-based environment before downloading the data. It’s important to note that limitations related to the size of the requested data can be a potential concern in this process. Alternatively, STAC simplifies the process, particularly when combined with cloud-ready formats. This lazily enables the creation of ESDCs. In the information-preserving data cube approach, platforms like ODC offer a comprehensive system for transforming original data into ESDCs and even provide mechanisms for storing the resulting ESDCs within the data cube infrastructureFootnote ¹¹. Noteworthy is openEO (Schramm et al., Reference Schramm, Pebesma, Milenković, Foresta, Dries, Jacob, Wagner, Mohr, Neteler and Kadunc2021), an API striving to connect multiple backends in a standardised way, including image collection providers (e.g. GEE) and information-preserving data cubes (e.g. ODC)Footnote ¹², to generate ESDCs.

3.1. Collection

Given that data providers frequently utilise diverse formats and protocols for data sharing, particularly in proportion to the multidimensional complexity of the data, the establishment of streamlined access mechanisms becomes imperative. Traditionally, File Transfer Protocol (FTP) servers have been used for data sharing. However, to enhance data discoverability and usability, data providers are increasingly adopting data stores that offer persistent and standardised data storage. Repositories play a vital role in this process by standardising metadata, enabling easy search and retrieval of assets through metadata queries. Recently, more and more data providers offer APIs to facilitate efficient querying of metadata and access to the data itself by adopting specifications such as STAC and enabling range requests for cloud-optimised data (e.g. Zenodo recently started to support HTTP range requestsFootnote ¹³).

The flexibility of these specifications enhances data interoperability by enabling the development of extensions that simplify data integration. For instance, the Electro-Optical STAC-extensionFootnote ¹⁴ has been created to facilitate the integration of multispectral remote sensing data by expanding the capabilities of STAC to accommodate specific requirements and metadata associated with this kind of data. Looking ahead, the advantages of data interoperability may potentially extend beyond the realm of raw source data, encompassing entire ESDCs. The datacube STAC-extensionFootnote ¹⁵ (currently in a “candidate” Extension Maturity level) has been developed with the primary objective of advancing the integration and interoperability of structured data representations like ESDCs within the STAC ecosystem. This effort aims to broaden the scope of opportunities for reusing ESDCs in new data processing pipelines.

Additionally, the efficiency of data access and collection is contingent upon data formats. GeoTIFF is arguably the most used and renowned data format for georeferenced raster data. This format adds a standard specification for the TIFF format that describes the spatial properties of the raster. It is widely used for EO products such as Landsat imagery. The need to operate in cloud environments has driven the development of cloud-optimised geospatial data formats. Consequently, the GeoTIFF format has evolved to the Cloud-Optimised GeoTIFF (COG)Footnote ¹⁶ format, enhanced to function efficiently in cloud environments through HTTP range requests. COGs offer several advantages over traditional GeoTIFFs, including reduced latency in data retrieval, faster visualization of large datasets, and a tiled structure that enables parallel processing. The significance of this format is underscored by its recent approval as an Open Geospatial Consortium (OGC) standardFootnote ^{17 ,} Footnote ¹⁸.

When the dimensionality of the data increases, formats such as NetCDF or HDF5 are typically used to encapsulate data and coordinate values. Tiling and chunking allow efficient access to big data arrays for both data formats. However, these formats are not inherently optimised for cloud environments. The Zarr specificationFootnote ¹⁹ addresses this limitation and can be used directly in cloud environments, offering several advantages over NetCDF and HDF5. Zarr enables more efficient chunk access for parallel processing, provides better support for distributed computing, and offers improved read and write speeds, particularly in cloud storage systems. Moreover, Zarr’s flexible chunking scheme allows for optimised data access patterns, and its simpler metadata structure facilitates easier data discovery and management. Additionally, specifications such as geo-zarrFootnote ²⁰ and the xcube dataset conventionFootnote ²¹ have been introduced to further enhance data interoperability and compatibility within the context of Earth system data.

3.2. Curation

Effective data curation stands as a critical anchor in the preparation of data for subsequent spatio-temporal processes and analysis via ESDCs (Marujo et al., Reference Marujo, Ferreira, Queiroz, Costa, Arcanjo and Souza2022). The transformation of raw data into Analysis-Ready Data (ARD) has emerged as an essential prerequisite across multiple initiatives. ARD ensures that data are readily amenable to queries, analysis, and application development. Notable instances of these initiatives include DE AfricaFootnote ²², DE AustraliaFootnote ²³, and the Brazil Data Cube (Ferreira et al., Reference Ferreira, Queiroz, Vinhas, Marujo, Simoes, Picoli, Camara, Cartaxo, Gomes, Santos, Sanchez, Arcanjo, Fronza, Noronha, Costa, Zaglia, Zioti, Korting, Soares, Chaves and Fonseca2020; Marujo et al., Reference Marujo, Ferreira, Queiroz, Costa, Arcanjo and Souza2022), among others.

The Committee on Earth Observation Satellites (CEOS) has precisely defined ARD as “satellite data that have been processed to a minimum set of requirements and organized into a form that allows immediate analysis with a minimum of additional user effort and interoperability both through time and with other datasets”Footnote ²⁴. In this definition, ARD also exhibit interoperability both across time and with other datasets (refer to Siqueira et al., Reference Siqueira, Tadono, Rosenqvist, Lacey, Lewis, Thankappan, Szantoi, Goryl, Labahn, Ross, Hosford and Mecklenburg2019 for an overview of the CEOS ARD for Land initiative, CARD4L). CEOS has established a comprehensive set of Product Family Specifications (PFS) tailored to various data groups, including Surface Reflectance, Surface Temperature, Polarimetric Radar, and more. These specifications undergo rigorous peer review processes across multiple satellite platforms, such as Landsat and Sentinel collections, to obtain CEOS ARD approval. It’s worth noting that there are ongoing efforts to develop additional PFS, including Interferometric Radar and LiDAR Terrain and Canopy Height. It is also worth noting that the OGC has recently addressed the CEOS ARD concept by forming a new Standards Working Group (SWG) to define a generic multi-part standard specifying a set of minimum requirements for geospatial products to be considered ARDFootnote ²⁵.

It is important to recognise that achieving an ARD level can extend beyond minimum standard specifications. Obtaining ARD often involves crucial preprocessing and data curation tasks that are tailored to the unique requirements of the application domain. For instance, in the context of EO data, these tasks may encompass but are not limited to, cloud and cloud shadow masking (refer to Skakun et al., Reference Skakun, Wevers, Brockmann, Doxani, Aleksandrov, Batič, Frantz, Gascon, Gómez-Chova, Hagolle, López-Puigdollers, Louis, Lubej, Mateo-García, Osman, Peressutti, Pflug, Puc, Richter, Roger, Scaramuzza, Vermote, Vesel, Zupanc and Žust2022 for a comprehensive intercomparison exercise of multiple cloud and cloud shadow masking methods), snow masking (e.g. Richiardi et al., Reference Richiardi, Blonda, Rana, Santoro, Tarantino, Vicario and Adamo2021), and the correction of Bidirectional Reflectance Distribution Function (BRDF) effects to derive Nadir BRDF Adjusted Reflectance (NBAR) values (e.g. Roy et al., Reference Roy, Zhang, Ju, Gomez-Dans, Lewis, Schaaf, Sun, Li, Huang and Kovalskyy2016).

3.3. Cubing

The concept of ARD may exhibit some subjectivity depending on the specific application. This subjectivity pertains to the data that populates an ESDC. In contrast, ESDCs inherently represent straightforward yet robust analysis-ready integrated entities, capable of simplifying a broad spectrum of analytical tasks (Baumann et al., Reference Baumann, Misev, Merticariu and Huu2019). An ESDC filled with ARD is often called an Analysis-Ready Data Cube (ARDC), a concept widely employed in DeepESDL. To generate an ARDC, the critical step involves aligning data onto a unified grid. Domain experts predefine this grid, and all data sources must conform. Furthermore, the efficacy of the ARDC processing is significantly influenced by the implementation of an optimal chunking strategy for this grid. This strategy must be tailored to facilitate efficient data processing across diverse analytical scenarios. For instance, analyses focused on temporal dynamics benefit from chunking strategies that preserve the temporal dimension, whereas spatial analyses or cartographic visualisations are optimised by maintaining spatial dimensions within chunks. In scenarios requiring multi-temporal spatial analysis, a hybrid approach combining both temporal and spatial preservation in chunking can be advantageous.

When the grid moves in the spatio-temporal domain, the varying spatio-temporal resolutions and coverage among multiple data sources require selecting adequate methods to fit the data into the predefined grid. Datasets with varying spatial resolutions and coverage must be resampled onto a standard spatial grid. This process often requires modifying the data (Cracknell, Reference Cracknell1998). While non-destructive algorithms such as nearest neighbours can preserve data values (at the cost of duplicating or ignoring values), significant differences in spatial resolution often require transformation through (non-) linear resampling methods, such as cubic convolution or advanced fusion techniques (Nikolakopoulos, Reference Nikolakopoulos2008). Complex AI methods can be employed to perform spatial transformations while preserving the quality of the measured variable (e.g. multi-image super-resolution algorithms, Michel et al., Reference Michel, Vinasco-Salinas, Inglada and Hagolle2022; Razzak et al., Reference Razzak, Mateo-García, Lecuyer, Gómez-Chova, Gal and Kalaitzis2023). Another often overlooked issue arises when dealing with extensive variables. In such instances, it is crucial to ensure that, for example, mass balances are not distorted in the process of creating new products.

It is important to note that the application of resampling methods, particularly in the context of generating Global ESDCs covering the entire planet (e.g. Mahecha et al., Reference Mahecha, Gans, Brandt, Christiansen, Cornell, Fomferra, Kraemer, Peters, Bodesheim and Camps-Valls2020), may introduce geometric distortions. Projecting global datasets onto a plane can distort the data in terms of area, distances, and angles (Snyder and Voxland, Reference Snyder and Voxland1989), posing challenges for subsequent analysis (cf. Section 5.1). This can be alleviated by using a Discrete Global Grid System (DGGS, Kmoch et al., Reference Kmoch, Vasilyev, Virro and Uuemaa2022). This kind of grid system seeks to minimise distortions, harmonise cell sizes and maintain consistent distances from neighbours. Defining standards and solutions for efficient chunk storage, subsetting, and integration into the ESDC framework will be a challenging future task. Still, it could lead to significant improvements in both the performance and accuracy of spatial algorithms.

In the case of Regional ESDCs (e.g. Estupiñan Suarez et al., Reference Estupiñan Suarez, Gans, Brenning, Gutierrez-Velez, Londono, Pabon-Moreno, Poveda, Reichstein, Reu and Sierra2021), which may cover entire continents, oceans, or administrative regions at various hierarchical levels, selecting an appropriate CRS is crucial to ensure minimal geometric distortion. On local scales, Local ESDCs (also referred to as mini cubes, Requena-Mesa et al., Reference Requena-Mesa, Benson, Reichstein, Runge and Denzler2021) cover smaller areas of interest (e.g. Walther et al., Reference Walther, Besnard, Nelson, El-Madany, Migliavacca, Weber, Carvalhais, Ermida, Brümmer, Schrader, Prokushkin, Panov and Jung2022), ideally characterised by high spatial resolutions ranging from sub-meters to meters. Using local ESDCs together with a local CRS enables to minimise distortions.

When dealing with datasets characterised by varying temporal grids, even if they share the same date-time units, irregular temporal grids may emerge. These discrepancies can introduce temporal gaps within the time dimension. In cases where datasets exhibit varying date-time units, especially when working with datasets featuring finer date-time units (e.g. daily records), it becomes necessary to aggregate them to align with a predefined coarser temporal grid (e.g. monthly records). While this process is straightforward for regularly sampled data, it can pose challenges for EO data with long revisit periods (e.g. Landsat data). These challenges can potentially introduce uncertainties during aggregation. Substantial gaps in EO data can have a detrimental impact on the accuracy and representativeness of the aggregated results. This concern is further exacerbated when additional gaps arise due to data disturbances, such as cloud and shadow interference.

3.4. Harmonisation

Additional post-processing of data variables may be necessary to address Earth system challenges. This entails further data curation to obtain a fully gap-filled, harmonised product with evenly spaced time steps. In alignment with the naming conventions established for EO data cubes by Frantz, Reference Frantz2019, we refer to a thoroughly harmonised ESDC as a highly ARDC (hARDC).

Data harmonisation is crucial to ensure the consistency and compatibility of variables obtained or generated using different methodological or technical approaches (Wulder et al., Reference Wulder, Hilker, White, Coops, Masek, Pflugmacher and Crevier2015). When discrepancies exist between data measurement or production methods, it can introduce inconsistencies that hinder subsequent analyses involving the specific variables (Vogeler et al., Reference Vogeler, Braaten, Slesak and Falkowski2018). To address this, one approach is to create separate variables that represent the same measured quantity, highlighting the differences between them. However, to enhance spatio-temporal resolution and coverage, harmonisation of variables is often necessary (e.g. harmonising reflectance values from Sentinel-2 and Landsat, Claverie et al., Reference Claverie, Ju, Masek, Dungan, Vermote, Roger, Skakun and Justice2018; Marujo et al., Reference Marujo, Carlos, da Costa, de Souza Arcanjo, Fronza, Soares, de Queiroz and Ferreira2023).

This can be achieved through simple methods that involve sampling data from the same spatio-temporal index in both variables to establish a direct conversion model (e.g. using matched observations to match Landsat 8 and Sentinel-2, Shang and Zhu, Reference Shang and Zhu2019). Alternatively, more advanced AI models can harmonise data by incorporating one or more additional variables (e.g. creating a global product of OCO-2 Sun-Induced Fluorescence, SIF, Li and Xiao, Reference Li and Xiao2019). This may require the development of an entire AI pipeline to extend a variable with newly available data or reconstruct it, especially in cases where the variable was not previously measured (e.g. reconstructing SIF from TROPOMI, Chen et al., Reference Chen, Huang, Nie, Zhang, Wang, Chen and Chen2022). In this sense, data harmonisation also encompasses projecting data in simulated future scenarios (e.g. projecting vegetation dynamics for the rest of the century, Mahowald et al., Reference Mahowald, Lo, Zheng, Harrison, Funk, Lombardozzi and Goodale2016). In addition, it is crucial to incorporate uncertainty metrics to facilitate accurate and reliable future analysis using the harmonised data variables (cf. Section 4.3).

Additionally, to effectively use algorithms that incorporate temporal structures, such as Recurrent Neural Networks (RNNs, Sherstinsky, Reference Sherstinsky2020), a regularly spaced and gapless time dimension is usually required. Hence, data from an irregular time dimension should be aggregated or interpolated to fit into a regular temporal grid (e.g. gap-filling Landsat reflectances on a monthly basis, Moreno-Martínez et al., Reference Moreno-Martínez, Izquierdo-Verdiguier, Maneta, Camps-Valls, Robinson, Muñoz-Marí, Sedano, Clinton and Running2020). A suitable predefined temporal resolution must be selected, and data must be gap-filled. Various gap-filling techniques, ranging from simple linear interpolation to more complex AI-based modelling approaches, can be employed to address this (e.g. using Long Short-Term Memory networks, LSTMs, Ren et al., Reference Ren, Cromwell, Kravitz and Chen2022). The choice of the gap-filling method depends on factors such as the data’s nature, the desired accuracy level, and the specific requirements of the analysis or application.

3.5. Transformation

Expertly crafted higher-order features often prove highly relevant for addressing Earth system challenges. These new features span a spectrum, encompassing operations that range from simple transformations of the original variables to the creation of entirely novel features derived from advanced AI models. Examples of such features include the computation of spectral indices derived from reflectance bands (Montero et al., Reference Montero, Aybar, Mahecha, Martinuzzi, Söchting and Wieneke2023), the extraction of frequencies through time series decomposition (Mahecha et al., Reference Mahecha, Fürst, Gobron and Lange2010), the creation of spatio-temporal compositions (e.g. Griffiths et al., Reference Griffiths, van der Linden, Kuemmerle and Hostert2013), summarising high dimensional dynamics (e.g. Kraemer et al., Reference Kraemer, Camps-Valls, Reichstein and Mahecha2020), and outputs generated by AI models (e.g. Brown et al., Reference Brown, Brumby, Guzder-Williams, Birch, Hyde, Mazzariello, Czerwinski, Pasquarella, Haertel, Ilyushchenko, Schwehr, Weisse, Stolle, Hanson, Guinan, Moore and Tait2022). To illustrate, consider a study focusing on climate extremes like heatwaves and droughts’ impact on the terrestrial biosphere. In such cases, calculating anomalies for critical variables (e.g. air temperature and soil moisture as proxies for heatwaves and droughts, with Gross Primary Production as the target biosphere variable) is pivotal (see Figure 3). Creating these novel features introduces a new dimension to distinguish between variable values corresponding to raw data and those representing anomalies. In line with the naming conventions introduced by Frantz, Reference Frantz2019, we designate an ESDC with higher-order features as a hARDC Plus (hARDC+).

Figure 3.

Abstract representation illustrating the connection between three Earth system variables in a hARDC+ (from top to bottom: anomalies in air temperature, soil moisture, and gross primary production). The arrows illustrate the interactions that can be modelled, e.g., predictive modelling (top to bottom) or interpretation (bottom to top), depending on the use case of interest.

3.6. Reuse

ESDCs, after generation and analysis, can either evolve into dynamic versions through continuous updates or become static ESDCs, serving as input for the generation of new ESDCs. In the first scenario, establishing a Continuous Integration (CI) pipeline becomes essential for automating ESDC updates. This pipeline can be scheduled to align with the release of new dataset versions, ensuring the ESDC remains current. However, this approach may prove inefficient for EO products that are delivered regularly (e.g. Sentinel-2 or MODIS) and that may be constantly reprocessed by data providers, releasing new versions with updated processing pipelines. In this case, the update schedule should align with the specific needs of the ESDC (e.g. monthly, semi-annually, or annually). In the second scenario, an automatic update of dataset versions is also feasible, eliminating the necessity to extend the ESDC to the most recent date. In either scenario, it is crucial to implement clear reproducibility and traceability practices to ensure the data integrity of future ESDCs.

As highlighted in Section 3.1, standardisation is pivotal in promoting fluid data interoperability within this context. OGC has recognised the growing significance of data cube approaches for geospatial data. OGC recently established the GeoDataCubes SWG to create an API that facilitates interoperability among various solutionsFootnote ²⁶. This standard covers a broad scope, explicitly including API functionalities for access and processing, exchange format recommendations, profiles, and a metadata model.

Additionally, cloud technologies have ushered in the development of data cube services that abstract the underlying file structures and formats, replacing them with APIs offering diverse processing functionalities and promoting interoperability. For instance, platforms like Sentinel HubFootnote ²⁷ serve as sources for ESDC generation through tools like xcube. Moreover, the openEO platformFootnote ²⁸ aims to provide an API that enables connections from multiple clients to various cloud backends using a unified API (Schramm et al., Reference Schramm, Pebesma, Milenković, Foresta, Dries, Jacob, Wagner, Mohr, Neteler and Kadunc2021). Approaches like these allow the tailored specification of ESDCs on-demand, with server-side processing relieving requesters of the complexities of the generation task. However, this convenience often comes with a trade-off, as the processing engine’s code basis, the processing environment, and the input data are not known to requesters. Any modifications to these specifications can result in different outcomes for identical requests to the data cube API, hindering a streamlined update of dynamic ESDCs and a transparent basis for reusing static ESDCs.

In contrast, less convenient but more transparent approaches fully document the ESDC generation process through “recipes”. These recipes contain versioned source code used for input data processing. Examples include the Pangeo ForgeFootnote ²⁹ (Stern et al., Reference Stern, Abernathey, Hamman, Wegener, Lepore, Harkins and Merose2022) and DeepESDL recipesFootnote ³⁰. Recipes, coupled with versioned input data and fully specified processing environments, enable practical reproducibility of resulting ESDCs. This approach supports the seamless updating of dynamic ESDCs when new data becomes available and provides transparency for incorporating static ESDCs into new datasets.

Ongoing efforts to enhance data lineage and provenance transparency are integral to the Copernicus Data Space Ecosystem. The development of the “traceability” serviceFootnote ³¹, currently in progress, is designed to empower users to trace all modifications to the data from its origin to its delivery to the end user, ensuring greater transparency and accountability in the ESDC life cycle.

3.7. Metadata generation

Traceability and self-explanatory power are essential aspects alongside the data values themselves. When an ESDC is generated, end users may access its description through various sources, including documentation that adheres to best practices for open data publishing within the Earth sciences. Such practices are supported by data journals (e.g. ESSDFootnote ³²) and scientific associations (e.g. AGU Open ScienceFootnote ³³), provided that the data producers have furnished comprehensive documentation. However, the data must carry its own encapsulated description in the form of metadata, which typically comprises a set of attributes represented as key-value pairs. This ensures the data contain relevant information about their characteristics, facilitating understanding and utilisation.

Metadata generation should begin at the initial stage of data collection, encompassing crucial information such as data descriptors (e.g. name, units, measurement methods and equipment, resolution), data transformations (e.g. resampling or interpolation methods), metadata transformations (e.g. renaming procedures, conventions conversion), and responsible producers (e.g. creator entity, data provider). This metadata generation process should be consistently maintained throughout the entire ESDC life cycle, documenting each step undertaken to derive the final product (e.g. storing the process graphs from openEO when using this platformFootnote ³⁴). This ensures comprehensive self-contained documentation of the history and processing of the ESDC.

While flexibility exists in metadata management, conventions are crucial when dealing with Earth system data. The Climate and Forecast Metadata Conventions (CF Conventions, Hassell et al., Reference Hassell, Gregory, Blower, Lawrence and Taylor2017), for instance, represent a comprehensive set of standards specifically designed for Earth system data stored in formats such as NetCDF (although they can be readily applied to other formats like Zarr). These conventions facilitate the creation of clear and detailed descriptions of data variables and coordinate dimensions. Furthermore, software like xarray (Hoyer and Hamman, Reference Hoyer and Hamman2017) can parse CF Conventions and leverage them for different ESDC processesFootnote ³⁵. Compliance with CF Conventions not only simplifies data sharing but also promotes interoperability among various data sources, ensuring that ESDCs adhere to established standards.

4.1. Adding factual knowledge via PIML

A great addition to Machine Learning (ML) modelling is combining the pure data-driven approach with factual knowledge of the system under investigation (Karniadakis et al., Reference Karniadakis, Kevrekidis, Lu, Perdikaris, Wang and Yang2021). PIML leverages domain knowledge (typically mechanistic models or differential equations) and flexible data-driven ML methods (typically neural networks). Consequently, PIML models respect physical boundaries more faithfully while being flexible enough to approximate arbitrarily complex non-linear functions from data (cf. discussion and references in Reichstein et al., Reference Reichstein, Camps-Valls, Stevens, Jung, Denzler and Carvalhais2019). ESDCs provide a unique structure to access multiple Earth system data streams, and the equation-based model describes the underlying process. Thanks to this ready availability of data and equations, exploring PIML models using a wide array of baseline models would be far easier and faster. The equations detailing a given variable could be added to the cube as a sub-field of the variable of interest in the same way that space and time are. The eventual implementation should consider the multi-platform and multi-language nature of ESDCs. As illustrated above, this requires a unified and robust approach that suits multiple use cases.

4.2. Sampling for AI in a complex system

Sampling on ESDCs is essential for learning the concrete interactions of drivers, spatial conditions, timing, and other determinants of specific processes and their implications. This involves strategically selecting a manageable subset from the ESDC. This selection process is particularly important for ML algorithms, as they rely on these subsets to establish a foundational understanding of the process to be analysed (Atkinson et al., Reference Atkinson, Stein and Jeganathan2022; Nikparvar and Thill, Reference Nikparvar and Thill2021). Pseudo-random sampling facilitates a broad and diverse data selection, while regionalised sampling uses specific patterns within the data for a more targeted analysis. The latter proves particularly advantageous when the research goal is to comprehend specific phenomena.

Constructing representative samples in Earth system processes must ensure an unbiased representation of the target variable. The multidimensional nature of Earth system processes poses sampling challenges across multiple variables. Consider, for instance, a study that aims at understanding the effects of climate extremes on the terrestrial biosphere using AI (Sippel et al., Reference Sippel, Reichstein, Ma, Mahecha, Lange, Flach and Frank2018). We know that climate extremes such as heatwaves, droughts, extreme precipitation, flooding, etc., are typically associated with multiple variables (Flach et al., Reference Flach, Brenning, Gans, Reichstein, Sippel and Mahecha2021). Additionally, such events can co-occur in unfavourable sequences, i.e., compounding heatwaves, droughts, or floods following droughts (Zscheischler et al., Reference Zscheischler, Martius, Westra, Bevacqua, Raymond, Horton, van den Hurk, AghaKouchak, Jézéquel and Mahecha2020). To understand such circumstances, one should consider the full spatio-temporal extended in all relevant dimensions, including derived meta-variables that describe the characteristics of these events, such as timing, duration, extent, and intensity (Flach et al., Reference Flach, Gans, Brenning, Denzler, Reichstein, Rodner, Bathiany, Bodesheim, Guanche and Sippel2017). Often, additional factors gain significance. For example, ecosystem responses to extremes vary in space depending on ecosystem conditions (Mahecha et al., Reference Mahecha, Gans, Sippel, Donges, Kaminski, Metzger, Migliavacca, Papale, Rammig and Zscheischler2017), land-cover types (Flach et al., Reference Flach, Brenning, Gans, Reichstein, Sippel and Mahecha2021), and associated impacts, e.g., on the carbon cycle (Sippel et al., Reference Sippel, Reichstein, Ma, Mahecha, Lange, Flach and Frank2018). Building suitable AI models that predict such impacts requires including static data (e.g. vegetation type).

Yet, the critical question is then: how to obtain adequate and balanced training and validation data? Earth system processes often involve rare events of extreme conditions, which may occur sporadically over time and space. This rarity can lead to imbalanced datasets, where certain classes of the target variable or ranges of continuous values are underrepresented. This also applies to ranges of continuous values in an imbalanced distribution. Imbalanced datasets affect the performance and generalisation of models trained on these samples. Achieving spatio-temporal representativeness in this context can be challenging. To train ML algorithms for effective recognition and understanding of these events, it is crucial to include additional sampling within the specific domains where these events occur. For example, when constructing datasets for global flood (Li et al., Reference Li, Dang, Li and Zhang2023) or cloud detection (Aybar et al., Reference Aybar, Ysuhuaylas, Loja, Gonzales, Herrera, Bautista, Yali, Flores, Diaz, Cuenca, Espinoza, Prudencio, Llactayo, Montero, Sudmanns, Tiede, Mateo-García and Gómez-Chova2022), the methodology involves initiating automatic sampling that covers a broad spectrum of ecosystem conditions. Simultaneously, manually selected events are introduced. This approach ensures a balanced representation of different classes in the dataset, thereby enhancing the algorithm’s capability to accurately predict such events. Figure 4 showcases a potential workflow where event detection is performed based on global ESDCs, and samples for high-resolution ML are extracted based on a systematic sampling strategy (e.g. Ji et al., Reference Ji, Fincke, Benson, Camps-Valls, Fernandez-Torres, Gans, Kraemer, Martinuzzi, Montero, Mora, Pellicer-Valero, Robin, Soechting, Weynants and Mahecha2024). Here, analysing land cover purity is an option (a relatively homogeneous land cover dominated by a single vegetation type allows for easier comparisons and subsequent analyses), as well as incorporating mixed land covers (which introduces heterogeneity and interactions among land cover types), providing more comprehensive information for model training.

Figure 4.

Abstract representation illustrating the process of sampling high-resolution mini cubes for further analysis by considering vegetation land covers and extreme events detected via a global ESDC. Note that sample mini cubes are specified in the spatial and temporal ranges of the detected extreme events (also considering their occurrence).

Finally, the selection of samples with the necessary data dimensions must align with the chosen algorithm. For instance, tabular-based algorithms like tree-based methods require 2-dimensional batches (sample and variable), which are selected as individual points from the spatio-temporal domain. DL methods like Transformers (Vaswani et al., Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser and Polosukhin2017) or RNNs, e.g., LSTMs (Hochreiter and Schmidhuber, Reference Hochreiter and Schmidhuber1997), which consider sequence (or positional) dependencies, require 3-dimensional batches (e.g. sample, timestep, variable) and extract samples usually as subsets of time series from the spatial domain. Convolutional Neural Networks (CNNs, LeCun et al., Reference LeCun, Boser, Denker, Henderson, Howard, Hubbard and Jackel1989) may be used with 4-dimensional batches (e.g. sample, height, width, variable) by taking spatial subsets or grids from the temporal domain. DL methods accounting for both spatio-temporal dependencies, such as 3DCNNs or Convolutional LSTMs (ConvLSTMs, Shi et al., Reference Shi, Chen, Wang, Yeung, Wong and Woo2015), require 5-dimensional batches (e.g. sample, height, width, timestep, variable) and extract samples as subsets of ESDCs.

4.3. Quantifying uncertainties

Uncertainty quantification is crucial to Earth science, providing a comprehensive assessment of the reliability and confidence associated with scientific predictions, model simulations, and observational data. Capturing and modelling uncertainty is a complex task as it arises from various sources such as data limitations, model approximations, and the inherent complexity of Earth system dynamics.

Uncertainty can be broadly categorised into two types: epistemic uncertainty and aleatoric uncertainty (Kiureghian and Ditlevsen, Reference Kiureghian and Ditlevsen2009). Epistemic uncertainty refers to the model’s confidence in its predictions and is related to the choice of model parameters. Techniques such as Bayesian inference or Dropout can estimate epistemic uncertainty (Srivastava et al., Reference Srivastava, Hinton, Krizhevsky, Sutskever and Salakhutdinov2014; Gal and Ghahramani, Reference Gal and Ghahramani2016). Bayesian methods assign probability distributions to model parameters, directly quantifying uncertainty. In DL, dropout-based methods create model ensembles by randomly dropping out units during training, providing a measure of uncertainty based on the variability among the ensemble members. While these techniques may not completely capture the underlying uncertainty due to assumptions made during modelling or training, they are practical and can be employed to estimate uncertainty. These methods can be computationally demanding and time-consuming, mainly when applied to real-time applications. However, advancements in cloud platforms and the Monte Carlo (MC)-Dropout technique have enabled reliable uncertainty estimates, even when working with massive amounts of data (Martínez-Ferrer et al., Reference Martínez-Ferrer, Moreno-Martínez, Campos-Taberner, García-Haro, Muñoz-Marí, Running, Kimball, Clinton and Camps-Valls2022). On the other hand, aleatoric uncertainty is associated with the noise or variability present in the data (e.g. data affected by natural variability, measurement errors, or other sources of noise) and cannot be reduced. Instead, it can be identified and quantified as part of the uncertainty characterisation.

ESDCs involving measurements or modelled data can be accompanied by associated uncertainty values. Data assimilation techniques are key in incorporating data into ESDCs while considering the associated uncertainties. Approaches such as Kalman filtering, variational data assimilation, or ensemble-based assimilation can effectively merge different data sources and quantify the resulting uncertainties (Mathieu and O’Neill, Reference Mathieu and O’Neill2008).

5.1. Geometric challenge on planet Earth

Most ESDCs covering the whole globe use a simple longitude-latitude plate-carrée projection, which fits the ESDC model very well. The approach also allows for efficient storage and subsetting of cubes to user-generated subsets corresponding to a bounding box. However, for advanced data analysis, equirectangular projections have two main drawbacks: 1) grid cells differing in latitude do not have equal area, and 2) the distances to nearest neighbours are not constant.

The first drawback introduces a sampling bias towards high latitudes in the data. This bias can affect the representativeness and accuracy of analyses (cf. Section 5.2), particularly for regions located closer to the equator. The most trivial cases are computations of scalars, like global means (e.g. Figure 5), which need to be weighted or approaches like principal component analyses that require area-weighed covariance matrices. Effects of this kind have been known for decades and are considered climate textbook knowledge (Storch et al., Reference Storch, Zwiers and Livezey2000). However, they remain a challenge, as we find them often ignored in ESDC analytics. Issues of this kind can be alleviated using area-weighted statistics, suitable for most linear algorithms, or by performing weighted sampling from grid cells. For advanced, often non-linear data science methods, considering the spherical geometry is much more challenging, and careful consideration is advised before naive applications are performed. Even when applying area-weighted statistics correctly, oversampled areas lead to unnecessary increases in storage requirements and computation time.

Figure 5.

Comparison of air temperature at 2 m from ERA5 with and without weighting on the global mean time series computation. This rather trivial example shows how radically wrong any computation can be if the spherical nature of planet Earth is ignored.

The second drawback is particularly significant when applying spatial convolutions or moving window operations. To address this, several approaches can be employed. One option is to use Spherical Harmonics for simple convolutions, providing a transformation that respects the spherical nature of the data (Wieczorek and Meschede, Reference Wieczorek and Meschede2018). Spherical Harmonics can also be used as coordinate embeddings for neural networks (Rußwurm et al., Reference Rußwurm, Klemmer, Rolf, Zbinden and Tuia2023). Another approach involves graph convolutions that consider varying distances to neighbours.

5.2. Spatio-temporal representativeness for an accurate model evaluation

Diagnostics on predictive modelling with ESDCs can be challenged by the representativeness and spatio-temporal structure of training data (Tobler, Reference Tobler1970; Meyer and Pebesma, Reference Meyer and Pebesma2021; Ploton et al., Reference Ploton, Mortier, Réjou-Méchain, Barbier, Picard, Rossi, Dormann, Cornu, Viennois and Bayol2020; Kattenborn et al., Reference Kattenborn, Schiefer, Frey, Feilhauer, Mahecha and Dormann2022). Assessing the accuracy of a prediction is statistically straightforward as long as reference data is available for the entire population or if a respective sample represents the spatio-temporal structure of the population (Wadoux et al., Reference Wadoux, Heuvelink, De Bruin and Brus2021; Brus, Reference Brus2021). However, many modelling tasks build on observations not representative of underlying temporal dynamics or an entire land surface variability (e.g. upscaling functional ecosystem properties from sparse and clustered FLUXNET sites). Such an imbalance in reference data may not necessarily lead to a bias in model coefficients (Pabon-Moreno et al., Reference Pabon-Moreno, Migliavacca, Reichstein and Mahecha2022). However, it may lead to inflated prediction accuracy estimates, given the commonly limited capacities of ML to extrapolate into the unknown, where the predictor-response relationship may depart (Ludwig et al., Reference Ludwig, Moreno-Martinez, Hölzel, Pebesma and Meyer2023). Thus, the accuracy assessment of a prediction estimated from clustered samples will not represent the factual accuracy of predictions beyond the reference data availability. This is critical for assessing the quality of a prediction itself and potential error propagation in subsequent analysis (Yates et al., Reference Yates, Bouchet, Caley, Mengersen, Randin, Parnell, Fielding, Bamford, Ban and Barbosa2018; Meyer and Pebesma, Reference Meyer and Pebesma2021; Mila et al., Reference Mila, Mateu, Pebesma and Meyer2022). It is advised that predictions should inform on the area of applicability (Meyer and Pebesma, Reference Meyer and Pebesma2021), i.e., the area in which the predictor-space is covered by the reference data and obtained predictive accuracies thereof are assumed to hold.

However, assessing the predictive performance of a model inside the area of applicability may be challenged by the spatio-temporal structure of the training and test data. Commonly, adjacent observations (both in time and space) are more similar (autocorrelated in space and time), and therefore accuracies determined from test observations near the training data will be more accurate (Roberts et al., Reference Roberts, Bahn, Ciuti, Boyce, Elith, Guillera-Arroita, Hauenstein, Lahoz-Monfort, Schröder and Thuiller2017; Dormann et al., Reference Dormann, M. McPherson, B. Araújo, Bivand, Bolliger, Carl, G. Davies, Hirzel, Jetz and Daniel Kissling2007). For instance, seasonal effects can inflate model performance when using test observations near training data in the temporal dimension. Dependence among training and reference data results in any case on optimistic estimates of model performance, meaning that such accuracies do not reflect the actual transferability of the model to unseen areas or time steps (Roberts et al., Reference Roberts, Bahn, Ciuti, Boyce, Elith, Guillera-Arroita, Hauenstein, Lahoz-Monfort, Schröder and Thuiller2017). For instance, Ploton et al. (Reference Ploton, Mortier, Réjou-Méchain, Barbier, Picard, Rossi, Dormann, Cornu, Viennois and Bayol2020) showed that ML-based models found accurate in the presence of spatial dependent training and validation data may learn spatial data structures instead of transferable relationships between a response (biomass) and the predictors (environmental variables and optical reflectance). This may not only lead to erroneous model transferability and extrapolation to new spatial or temporal domains but also prevent an adequate interpretation of model functioning and attribution to variables and processes (Sweet et al., Reference Sweet, Müller, Anand and Zscheischler2023). Therefore, model performance and interpretation should be performed by minimising spatio-temporal dependence of observations via cross-validation strategies (cf. Roberts et al., Reference Roberts, Bahn, Ciuti, Boyce, Elith, Guillera-Arroita, Hauenstein, Lahoz-Monfort, Schröder and Thuiller2017; Meyer et al., Reference Meyer, Reudenbach, Hengl, Katurji and Nauss2018; Ploton et al., Reference Ploton, Mortier, Réjou-Méchain, Barbier, Picard, Rossi, Dormann, Cornu, Viennois and Bayol2020; Kattenborn et al., Reference Kattenborn, Schiefer, Frey, Feilhauer, Mahecha and Dormann2022).

6.1. Computing resources

The data size and available computing resources determine data processing feasibility throughout the ESDC life cycle. Computing resources vary from a single laptop to a local cluster with multi-threaded or distributed processing capabilities and can extend to cloud computing environments composed of multiple clusters. Modern computers are equipped with high-speed Solid-State Drives (SSDs) featuring fast random access and the potential for multiple Gigabytes per second throughput. However, the challenge lies in their limited capacity. In data centres, this is solved by using arrays of disks, but this introduces additional challenges, including latency, throughput, reliability, and security. Computation on local systems typically involves single-threaded or lightly multi-threaded computations with a higher level of interactivity. In High-Performance Computing (HPC) environments, the software operates in a multi-threaded or multi-core manner and is usually installed by a local system administrator. HPC environments are well-suited for extensive processing tasks but offer reduced interactivity due to the involvement of job schedulers for managing computation resources. Cloud computing environments offer a promising solution for managing vast amounts of Earth system data. These environments can be further improved in terms of scalability by utilising technologies like Kubernetes and Argo, which allow for specialised workflows. Platforms such as GEE, the European Open Science Cloud (EOSC)Footnote ³⁶, Google ColaboratoryFootnote ³⁷, Amazon SageMakerFootnote ³⁸, DeepESDLFootnote ³⁹, Copernicus Data Space Ecosystem (CDSE)Footnote ⁴⁰, and KaggleFootnote ⁴¹ provide opportunities for efficient data storage, processing, and collaboration in scientific research. However, it is essential to note that these platforms often impose certain limitations on the users. These limitations include storage capacity, computational resources, available tools for ESDC management, access permissions, and usage restrictions.

6.2. Software capabilities

In the context of managing ESDCs, diverse tools are available. Here, we present a compendium of useful tools for processing Earth system data within the ESDC life cycle in three prominent programming languages: Python, R, and Julia.

Python, arguably the most used language for ESDC management, offers xarray with labelled multidimensional arrays (Hoyer and Hamman, Reference Hoyer and Hamman2017), built on top of numpy (Harris et al., Reference Harris, Millman, van der Walt, Gommers, Virtanen, Cournapeau, Wieser, Taylor, Berg, Smith, Kern, Picus, Hoyer, van Kerkwijk, Brett, Haldane, del Ro, Wiebe, Peterson, Gérard-Marchant, Sheppard, Reddy, Weckesser, Abbasi, Gohlke and Oliphant2020), and supporting on-disk reading and parallel processing via dask (Rocklin, Reference Rocklin, Huff and Bergstra2015) (a Python library for parallel computing, enhancing array objects by employing data partitioning into chunks and employing dynamic task scheduling). Multiple tools are tailored to construct and process xarray datasets, which represent ESDCs. For data collection, rasterio (Gillies et al., Reference Gillies2013), rioxarray Footnote ⁴², satpy (Raspaud et al., Reference Raspaud, Hoese, Lahtinen, Holl, Finkensieper, Proud, Dybbroe, Meraner and Strandgren2023), or EOreader (Maxant et al., Reference Maxant, Braun, Caspard and Clandillon2022) are instrumental for reading GeoTIFFs and COGs, returning xarray objects. Xarray excels in reading NetCDF files and cloud-based data via zarr as dask-arrays. Vector data can be converted into xarray objects using geocube (Snow et al., Reference Snow, Taves, SlapDrone, Abdalla, Pierrick and Bell2023). Data sourced from STAC catalogues can be sought through pystac-client and directly transformed into xarray objects via stackstac Footnote ⁴³, odc-stac Footnote ⁴⁴, or cubo (Montero et al., Reference Montero, Aybar, Ji, Kraemer, Söchting, Teber and Mahecha2024). These tools support data collection and immediate cubing, including the temporal dimension. GEE enables data retrieval as numpy arrays through its API, which can be directly converted into xarray objects using Xee or wxee. GEE’s API (Gorelick et al., Reference Gorelick, Hancher, Dixon, Ilyushchenko, Thau and Moore2017) and extensions (Montero, Reference Montero2021) allow data curation before cubing. Xcube has various data stores for data acquisition and xarray object generationFootnote ⁴⁵. XDGGS (Kmoch et al., Reference Kmoch, Bovy, Magin, Abernathey, Coca-Castro, Strobl, Fouilloux, Loos, Uuemaa, Chan, Delouis and Odaka2024) simplifies working with different DGGS in xarray. The curation, harmonisation, and transformation stages, being subjective and application-dependent, can be accomplished through xarray or numpy processing. Libraries like scipy (Virtanen et al., Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng, Moore, VanderPlas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris, Archibald, Ribeiro, Pedregosa and van Mulbregt2020), built on top of numpy, offer additional resources leveraging ESDCs as multidimensional arrays. The analysis phase leverages a plethora of tools. ESDCs as multidimensional arrays are compatible with numpy, scipy, and related tools. Moreover, ESDCs represented as tensors interface effectively with tensorflow (Abadi et al., Reference Abadi, Agarwal, Barham, Brevdo, Chen, Citro, Corrado, Davis, Dean, Devin, Ghemawat, Goodfellow, Harp, Irving, Isard, Jia, Jozefowicz, Kaiser, Kudlur, Levenberg, Mane, Monga, Moore, Murray, Olah, Schuster, Shlens, Steiner, Sutskever, Talwar, Tucker, Vanhoucke, Vasudevan, Viegas, Vinyals, Warden, Wattenberg, Wicke, Yu and Zheng2016) or pytorch (Paszke et al., Reference Paszke, Gross, Massa, Lerer, Bradbury, Chanan, Killeen, Lin, Gimelshein, Antiga, Desmaison, Köpf, Yang, DeVito, Raison, Tejani, Chilamkurthy, Steiner, Fang, Bai and Chintala2019). Furthermore, developments that aren’t designed for direct ESDC use can also be leveraged using tensors in the representation of ESDCs (e.g. torchgeo, Stewart et al., Reference Stewart, Robinson, Corley, Ortiz, Lavista Ferres and Banerjee2022, GeoTorchAI, Chowdhury and Sarwat, Reference Chowdhury and Sarwat2022, pytorch-metric-learning, (Musgrave et al., Reference Musgrave, Belongie and Lim2020) and TorchIO, (Pérez-García et al., Reference Pérez-García, Sparks and Ourselin2021)).

R, a widely used programming language for statistical analysis, has assumed increasing significance in geospatial data processing and management. Raster data sourced from image collections can be managed seamlessly, progressing from data collection to study, with the assistance of libraries like raster Footnote ⁴⁶ or its more recent counterpart, terra Footnote ⁴⁷. ESDCs can be collected and analysed through dedicated tools like gdalcubes (Appel and Pebesma, Reference Appel and Pebesma2019) and stars (Pebesma and Bivand, Reference Pebesma and Bivand2023). Regionalised sampling using geospatial data can be conducted using stpp (Gabriel et al., Reference Gabriel, Rowlingson and Diggle2013) and spatstat (Baddeley et al., Reference Baddeley, Rubak and Turner2015). Recent developments have introduced the capability for lazy on-disk reading of Zarr filesFootnote ⁴⁸. Furthermore, data can be sourced and cubed directly from STAC catalogues using rstac (Simoes et al., Reference Simoes, de Souza, Zaglia, de Queiroz, dos Santos and Ferreira2021b) in combination with gdalcubes. Another comprehensive package for ESDC management is sits (Simoes et al., Reference Simoes, Camara, Queiroz, Souza, Andrade, Santos, Carvalho and Ferreira2021a), offering an end-to-end solution that additionally includes various tools for AI-related tasks, encompassing sampling, tuning, prediction, and the computation of uncertainty values.

Julia, a high-speed programming language, has gained popularity in scientific computing, making it an excellent choice for processing the large volumes of data found in ESDCs. Julia offers tools that cover crucial parts of the ESDC life cycle. These tools include YAXArrays.jl Footnote ⁴⁹ and Rasters.jl Footnote ⁵⁰ for multidimensional labelled array operations, GriddingMachine.jl (Wang et al., Reference Wang, Köhler, Braghiere, Longo, Doughty, Bloom and Frankenberg2022) for data acquisition, and experimental libraries like STAC.jl Footnote ⁵¹ for data discovery within STAC catalogues. For analysis, Julia provides specialised tools such as EarthDataLab.jl Footnote ⁵² for the direct processing of the Earth System Data Cube (Mahecha et al., Reference Mahecha, Gans, Brandt, Christiansen, Cornell, Fomferra, Kraemer, Peters, Bodesheim and Camps-Valls2020). Moreover, data distortions introduced during the cubing process can be addressed using libraries like OnlineStats.jl Footnote ⁵³ (Day and Zhou, Reference Day and Zhou2020) and WeightedOnlineStats.jl Footnote ⁵⁴ (Kraemer et al., Reference Kraemer, Camps-Valls, Reichstein and Mahecha2020). Julia’s ecosystem also includes ML tools like Flux.jl (Innes, Reference Innes2018), DiffEqFlux.jl (Rackauckas et al., Reference Rackauckas, Innes, Ma, Bettencourt, White and Dixit2019), and ReservoirComputing.jl (Martinuzzi et al., Reference Martinuzzi, Rackauckas, Abdelrehim, Mahecha and Mora2022), enabling advanced data analysis, including novel techniques like PIML.

6.3. Scalability obstacles

The size of ESDCs poses several challenges for analysis. Generally, in most programming languages for data science (e.g. Python, Julia, R), data has to be completely loaded into memory before calculating a simple statistic (e.g. median). However, ESDCs often surpass the memory limit, hindering computations or resulting in significant slowdowns due to frequent disk read-write operations. Instead, users can apply specialised algorithms that calculate statistics iteratively (Welford, Reference Welford1962; Schubert and Gertz, Reference Schubert and Gertz2018). memory algorithms allow the user to track statistics (e.g. mean, sums, and standard deviations) iteratively. They give the user complete control (and responsibility) over the order of the data reads. Because of the spherical nature of the Earth and the resulting differences in the area covered by pixels, these computations require weighted versions of the statistics (cf. Section 5.1). Errors arising from floating-point arithmetic must be minimised, including the potential for catastrophic cancellation (Kahan, Reference Kahan1965; Goldberg, Reference Goldberg1991).

Often, analyses can be performed independently on timesteps, maps, or any other discrete chunks of an ESDC (e.g. dimensions, periods, spatial slices). First, users split the data into those chunks, and then apply the transformation. In the end, users combine the elements back together into a new ESDC (see Figure 6). Many analyses can be expressed in terms of split-apply-combine (Wickham, Reference Wickham2011; Mahecha et al., Reference Mahecha, Gans, Brandt, Christiansen, Cornell, Fomferra, Kraemer, Peters, Bodesheim and Camps-Valls2020), such as calculating mean seasonal cycle maps from a time axis to a day-of-year axis, or a global mean temperature time series that collapses latitude and longitude into a scalar value per timestep. This method is also known as map-reduce in distributed data processing. Still, in contrast, it is made for array-like or tabular data (and the reduce step always consists in concatenating the results of the map step, cf. Wickham, Reference Wickham2011). Implementations of split-apply-combine can trade-off between memory consumption and performance by adjusting the amount of data being loaded into memory simultaneously. They may also take advantage of parallel reading, processing, and writing of data, which is especially important if the data is not stored on local storage but on object stores with high access latency.

Figure 6.

Split-apply-combine: split an ESDC along arbitrary axes, apply a function $ f $ to each sub-cube, and then combine the results along the same axes that have been used to split the original ESDC.

Storage in the form of compressed chunks typically employed by ESDCs, where reading a single element requires loading an entire chunk into memory, presents an opportunity for optimising sampling during ML training. Reading points individually is inefficient, as sampling two points from the same chunk necessitates reading the entire chunk twice. To mitigate this, reordering the points within a batch enables reading points from the same chunk jointly, reducing the number of reading operations. Adopting this approach makes it possible to limit the need to read the entire ESDC only once per batch, optimising the data access process.

Ensuring that scalability obstacles are transparent for end users during Earth system data analysis is essential. While experienced users may be able to address scalability issues effectively, less experienced users may struggle with the process if it is not fully transparent. It is important to provide a user-friendly interface that hides the complexities of scalability, allowing users to focus on their analysis tasks. Not all users can access sufficient computing resources for scaling processes, resulting in additional processing costs. Therefore, providing accessible and cost-effective solutions for scalability, such as cloud-based platforms, is crucial to enable a broader range of users to harness the benefits of scaling in Earth system data analysis.