B.6.1.1.  Pixel-in-X

A recurring discussion on geo raster data is whether the pixels are assumed to sit in the center or a corner of the cell. The root cause is a gross misinterpretation of diagrams.

In a 2D diagram the coordinate axes determine the position of points drawn in the diagram. Auxiliary lines go through the integer values or some other meaningful line spacing. Crossing auxiliary lines encloses areas, but in mathematics this has no meaning. If points are to be plotted at coordinates where auxiliary lines cross, then the points will invariably be plotted at the cross-section of two lines.

In programming, arrays (which are used to ultimately store the pixels) are diagrammatically usually shown as a lineup of boxes where each box symbolizes a storage area that holds some value. The borders in this diagram have no meaning.

Figure B.3 puts both types of diagrams next to each other. Note how in the left Cartesian diagram numbers (i.e., coordinates) are associated with the lines whereas in the right array diagram numbers (i.e., storage addresses) are associated with the boxes. Combining both concepts results in the notion of pixels that (i) have been acquired for some geographic point (left interpretation) and (ii) have been stored in some memory cell.

This has caused confusion among some geo scientists who started mixing both concepts and coming up with the idea that a pixel could sit in the “center” of the cell, leading (with regular grids) to a half-pixel offset of the pixel’s direct position. Consequently, the scientists differentiate between “pixel-in-center” and “pixel-in-corner” often making some silent assumption that the corner under consideration is the upper-left one while for d dimensional data 2d corners exist.

This interpretation is factually wrong. On the level of abstraction of gridded geographic data there is no concept of a “cell” in the array sense, just direct positions (in coverage terminology) sitting at the given real-world coordinates. Each sensor delivers data for the given coordinates, meaning at the crossing of the auxiliary lines for these coordinates. If the value was assumed to sit, say, “between x_i and x_i+1”, then it would have the direct position coordinates given by (x_i + x_i+1)/2. On a side note, some formats offer slots to store a local reference frame which could capture such an offset.

Figure B.3 — Cartesian diagrams (left) versus array symbolization (right)

The situation is aggravated by the lack of indicating the pixel-in-X assumption in some data formats. Combining datasets (including a background map) with different assumptions can lead to unwanted half-pixel shifts, potentially leading to dramatic deviations in an analysis. Discussions of how to handle this appear periodically on the Web, such as on gdal-dev where a contributor recently looked at the NetCDF format and Proba-V satellite data observing that wrong interpretation “could be dangerous or at least inconvenient. In my case I have a Python code that does zonal statistic and when a geometry is converted to a little raster binary mask, it does not match the same area in the ProbaV data and the results is a crash in Fiji because ReadRaster goes overboard or lost of NDVI data in other regions. The half pixel shift grows as you go east.”

This offset issue becomes even more spectacular when applied beyond the normally considered horizontal axes. In x/y/z coverages including a height or bathymetry axis there is no generally accepted convention that would determine a half-voxel offset upward or downward.

Similar questions arise on the temporal axis. Consider an example where some coverage with temporal resolution 1 year provides weather observations every January. Naturally, it would be expected that data “sit” in January. A pixel-in-center assumption would have a half-resolution offset, so 6 months. Hence, the weather observations would be attributed to June, rather than January.

Fortunately, the pixel-in-X discussion has not yet reached vertical and temporal axes and remains an item of discussion and confusion on the horizontal axes. Another difficulty is when non-regular grids come into play. Gridded coverages can be of the regular structure assumed so far (Figure B.4 left), or can be irregular of several types (Figure B.4 center and right).

Figure B.4 — Sample regular and irregular grid coverage types [26]

The recommendation is to improve education on this misconception and always assume a pixel sitting at the coordinate indicated, with no offset to anywhere. Course material should be developed and distributed widely which discusses pixel-in-X in particular on 4D spatiotemporal coverages.

B.6.1.2.  Pixel-is-X

Another issue centers around the question whether a pixel is considered a point (named “pixel-is-point”) or an area “around” that point (named pixel-is-area”). Typically, in regular grids this area would stretch halfway to the neighbor positions (Figure B.5). Motivation is that optical sensors average light values over an area while elevation data represent a sample valid at exactly (and only) at the measured point.

Figure B.5 — Pixel-is-area perception

There is no commonly defined handling. Users are expected to understand and interpret data appropriately. For example, the United States Geological Survey (USGS) points out [34] that Landsat uses a pixel-is-point interpretation and warns that some tools assume pixel-is-area designation. According to the above argument, however, Landsat should have a designation pixel-is-area.

To investigate this further, the characteristics of a sensor are considered. Details like viewing angles of a sensor are ignored, assuming an orthorectified, radiometrically corrected situation for simplicity. Likewise, the sensor instrument structure and its common variants such as pushbroom, whiskbroom, array, etc., are ignored.

Regardless of whether sensors are single or combined (e.g., CCD arrays), each single sensor has a so-called field of view (FOV) which is defined as the maximum angle of view that a sensor can effectively collect signals (such as photons). The width on the ground corresponding to the FOV is called the swath width.

A sensor will collect signals over some time over its FOV, resulting in the final pixel value. This FOV, due to the physical nature of the sensor, will be a cone which is narrow close to the sensor and wider as distance increases. The area from which a sensor captures signals can be described through an integral over the FOV area, which is circular or more generally elliptic.

In any case, the FOV cone does not “know” about the neighbor sensors, and by no means can it be assumed that signal capture stops midway between the FOV centers (Figure B.6 left), or conversely that every point in the swath width contributes to some pixel in the sensor array (Figure B.6 right). Even in the case of no overlapping and maximum FOV relative to sensor distance, it noticed that there are areas that do not contribute to any pixel.

Figure B.6 — Various FOV situations

Therefore, the assumption is that a pixel can always represent an area of signal collection (pixel-is-area), even if the area shrinks to a singular point (pixel-is-point). The main point is that this differentiation should not lead to a substantially different treatment, such as introducing a map offset — at least not on the level of abstraction discussed and notwithstanding particular sensor characteristics to be captured.

The main difference remaining is interpolation: The DEM should not be interpolated while it is perfectly fine for the optical sensor. In this opinion, it is captured already in the interpolation information of the coverage where the DEM would have no interpolation associated whereas the optical sensor is fine with various sorts of interpolation.

At a next level of differentiation, interpolation might be restricted to the FOVs of the sensors, but that information is not captured in the current standards and would possibly lead to a massive increase in complexity.

Finally, an observation is that pixel-is-x exclusively focuses on horizontal space. No guidance was found on how to apply these concepts on vertical and temporal coordinates. Additionally, reservations from above concerning shifting by a half-resolution offset apply.

Notably, in the OGC there is some controversy on this topic[35].

Further investigation into the technical details of where in the coverage extent interpolation should be allowed is recommended which should not stop at area versus point discussions but also drill more deeply into theory and practice. In any case, it is strongly recommended to abandon the idea of half-pixel shifts.

B.6.1.3.  Units of Measure

A long-standing, yet not satisfactorily solved issue is the unit of measure (uom). CIS tentatively does not make any coverage-specific assumption about its representation. After all, uom is an overarching concept which should not be made specific to coverages. A string-valued attribute for identification, as adopted from OGC Sensor Web Enablement (SWE) Common where attributes code and href, is suggested. SWE requires either a code attribute using UCUM or a reference to an external unit definition.

However, diving into the corresponding XML schema, it appears that the underlying uom code type, UnitReference, is not defined. Instead, there is a UnitReferencePropertyType type whose attribute code is of type UomSymbol. Even in the examples provided with the SWE Common Standard, a wide range of usage is found:

• use of the code attribute for a symbol, name, abbreviation thereof, or a formula: code=”%”, code=”mbar”, code=”deg”, code=”km/h”; and

• use of the xlink:href attribute for a reference to a concept or an abbreviated unit name (for which, however, no registry exists).

xlink:href=”http://www.opengis.net/def/uom/ISO-8601/0/Gregorian”
xlink:href="Cel"

In CIS, the range type — which is adopted from SWE Common — contains a corresponding element which is supposed to contain a UCUM unit or a URI, as per OGC convention. A URL would act as an identifier of some unit, and with a suitable mechanism behind this would allow even automatic transformation between units, say between feet and meters. Similarly, in the domain definition CIS remains agnostic and just requires units — such as geographic degrees or km, or temporal timestamps — which are understood in the context used.

Actually two issues can be identified. First, an incoherence in names can be observed. This happens when abbreviations are used. Even on very basic geographic units, uniformity is lacking as the following two random examples demonstrate.

• OGP at some time changed “long” for longitude to “lon” in the EPSG axis abbreviation to have a three character abbreviation like “lat”. This led to major issues in the handling of existing and new data as software usually did not automatically recognize that “long” and “lon” both denote the same axis. In the end, the solution in OGC was to allow in a coverage to freely name the axis labels, not tied to the name used in the EPSG CRS definition. These axes could be identified by position in the axis sequence — for example, CRS axes “Lat Lon h ansi” in a coverage can be named “lat long height date”. This decoupling allowed tools to keep legacy “long” and also to liberate applications from using cumbersome OGC names such as “ansi” for dates.

• Another issue of the same kind, which remains unsolved today, is the use of both “deg” and “degree” for the unit of measure along horizontal geographic axes.

Second, coverages are lacking a framework which would allow automatic conversion across the manifold units in practical use worldwide.

Obviously, there is a lack of standardization for such information, leading to nonhomogeneous or missing information in operational data. This has been observed already earlier [36]. Relevant candidates are inspected.

UCUM (Unified Code for Units of Measure) is a code system intended to include all units of measures being contemporarily used in international science, engineering, and business [37]. Units can be constructed following a mini language. For example, radiation per area might be expressed as “W/cm2” for Watts per square centimeter. If no unit is to be provided, by CIS convention, uom is “1” which in UCUM is expressed as “10^0”. Notably, for complex physical phenomena UCUM can become quite involved, as the radiance definition “W.m-2.sr-1.nm-1” in one CIS schema example shows [45].

The European Environmental Agency (EEA) has published a dictionary with concepts relevant for environmental monitoring and land management [38]. This hierarchical classification system uses the Simple Knowledge Organization System (SKOS), a common data model for sharing and linking knowledge organization systems via the Web [39]. According to the website, this RDF-based data dictionary was last updated in 2015, so does not appear well maintained. Further, the power of ontologies seems not fully utilized as the main connectivity between terms is the hierarchy, without a fully-fledged explicit ontology network.

RAINBOW is the OGC Definitions Server [40] which offers information about concepts relevant in the standards ecosystem, with units of measure, CRSs, etc. While the service offers a wide range (including, for example, planetary science CRSs) RAINBOW is relatively new. For example, for temperature there is Kelvin, but neither Celsius nor Fahrenheit are available, and the EPSG section is completely empty. UCUM is referenced, but the corresponding concept is not connected to anything. Most important for this discussion, no information is provided that could be used for some automatic conversion, such as from seconds to minutes or meters to feet. A discussion has been launched recently, but subsequently has been closed again [41].

An alternative is QUDT [46], a unified architecture for the conceptual representation of quantities, quantity kinds, units, dimensions, and data types. It appears that QUDT is substantially more comprehensive than the uom ontology started by OGC. For example, temperature is linked to Kelvin, Fahrenheit, and Celsius, among others. DEG_C (Degree Celsius) contains information on the derivation from the canonical (Kelvin) unit indicating a conversion multiplier of 1.0 and a conversion offset of 273.15, in a machine-readable format. This indeed allows for an automatic conversion.

A potential downside of QUDT is the lengthy URIs which are unwieldy, at least for human users (and even more so with composite CRSs, see discussion of multi-dimensional CRSs). One could argue, though, that only the developers see such URIs. In operational applications, there is only machine-to-machine communication where the length of the notation does not matter.

Notably, QUDT is not a standard, but a research outcome maintained by a small (but rather active) group. Custom datatypes for measurements as well as on-the-fly support of arbitrary custom datatypes are implemented in an Apache Jena fork.

Recently, another RDF ecosystem of practically relevant datatypes has been suggested [47], including a microformat for UCUM syntax next to energy, force, pressure, speed, temperature, time, etc., which combines the exact notation of UCUM with the reasoning capabilities as in QUDT. An implementation based on the Apache Jena reasoning framework is available. In future this might open the door for an automatic conversion across atomic and composite units.

A further recent activity has been launched by of CODATA with its Digital Representation of Units of Measurement (DRUM) task group [48]. The Bureau International des Poids et Mesures (BIPM) [49], an international organization established by the Metre Convention and managing body of the International System of Units (SI) and the international reference time scale (UTC), has started, during Testbed 19, a Forum on Metrology and Digitalization with one aim being a SI Digital Framework. Both CODATA and BIPM goals appear similar to what is discussed in this ER Annex, but no results are available yet. At the University of North Florida, a prototypical Units of Measure Interoperability Service [50] has been established aiming at automatic units translation.

Comparing the approaches discussed so far reveals the following findings.

• Common sense names and abbreviations are not suitable for unit understanding and automatic conversion.

• UCUM is exact and can be parsed for automatic handling by introducing a microsyntax which is nicely compact but can become unwieldy for human users.

• OGC’s ontology approach is going the right direction but is not (yet) usable in its current preliminary state.

• QUDT is amenable to reasoning as well as to conversion arithmetics.

Bottom line, it is recommended adopting QUDT for use in the range type — best as a normative requirement, but at least as a good practice. OGC might be the right body for that. Also, OGC could host the specification and service as part of OGC’s “community standards” process which has become popular recently.

B.6.1.4.  Tiling

Partitioning of large arrays into smaller sub-arrays, also known as tiling or chunking, is a common technique for storing large arrays. Queries requiring some region of a gridded coverage will need to retrieve the required tiles, cut out the intersection with the region, and reassemble the parts into the result.

A comprehensive analysis of tiling in general has been provided by Furtado [51]. Particular patterns were categorized and cast into several parametrized strategies for easier handling by the tiling responsible, normally the data manager (Figure B.7). In rasdaman, these strategies have been implemented in a storage layout clause extending the insert and update statements. The engine will use the tiling directives for all new incoming data arranging the data according to the pattern specified. Via an update statement, a repartitioning is possible, something useful if new insights about the access patterns have emerged.

Figure B.7 — Sample tilings, after Furtado: from left, regular, aligned, non-aligned, area-of-interest strategy

The goal of tiling is always to minimize the number of disk accesses. Within limits it is of less importance how large the tiles are, but every I/O request costs. What is important is to minimize those I/O requests. If tiles are too small (like 4kB at some time used by Oracle and Esri, and about 100×100 pixels used in the 2D tiling of PostGIS Raster [52]) too many loads need to be performed, leading to significant performance losses as verified experimentally by Furtado. If tiles are too large, on the other hand, too much data are loaded and then dropped while cutting out the result. As such, tiling becomes a tuning parameter giving the administrator complete freedom to define any partitioning, or simply rely on some reasonable defaults.

Array databases hide the underlying storage structure of the datacubes, thus relieving users from a potentially challenging technicality which has proven useful in practice for simplifying access APIs and client codes. Some tools, though, require the user to reassemble the tiled and clipped data. In OGC standardization attempts exist, though, to make the tiling explicit.

OGC Web Map Tile Service (WMTS) “trades the flexibility of custom map rendering for the scalability possible by serving of static data (base maps) where the bounding box and scales have been constrained to discrete tiles”[19]. This reduced functionality is justified by easier implementation: “The fixed set of tiles allows for the implementation of a WMTS service using a web server that simply returns existing files.” This tiling, though, is offered only for 2D x/y, any other potential dimensions a datacube can have are not considered. Elevation is considered in the 3D Tiles standard [53], however not time. On a side note, geographical parameters also are hardwired: Angles are fixed to radians, distances to meters.

OGC Abstract Topic 22 [54] specifies models for 2D tilings which, without further justification, are said to extend into higher dimensions. The 3D Tiles standard [53] is yet another specification addressing tiling. Further, in ISO a coverage tiling standard proposal has been submitted recently.

Notably, Abstract Topic 22 emphasized that a tile internally obeys a “homogeneity constraint” of following a single tessellation rule. A tile set making up one particular object is assumed to have a common CRS, uom, origin, and extent for all its tiles. In contrast, in the coverage world the coverage is the unit of homogeneity. All these tiling concepts and APIs have in common the introduction of a concept at the user level which is not required by the scientific task (such as EO analysis) and substantially complicates handling coverages.

One argument is that visual map clients request data in a tiled way and, therefore, providing tiles is natural. However, storage tiling and delivery tiling are two separate concerns. A client might get data in tiles, but a spatiotemporal access API should always allow user-selected arbitrary bounding boxes and not constrain to some storage pattern the server finds useful. Finally, while visual clients like to access in tiles of their own chosen sizes, this certainly does not hold for analysis where, for example, the code for processing a 2D array should be a simple nested loop on a single, contiguous array. Everything else creates additional complexity, programming needs, and opportunities for introducing bugs — in short: it is not analysis-ready.

Notably, CIS supports tiled storage, but that is part of the internal data organization — neither WCS GetCoverage nor WCPS queries are aware of such internals. Both operate regardless of what the tiling of a coverage is like.

Therefore, it is ultimately claimed that for analysis-readiness tiling should be opaque and hidden from the user. If deemed necessary to standardize, then a separate tuning API might be defined which could fit well into the coverage administration API, WCS-T [55].

B.6.1.5.  Coverage Metadata

The coverage metadata element contains any ancillary information beyond the canonical metadata fixed in the coverage structure. The corresponding XML schema part of CIS is as follows:

<complexType name="MetadataType">
    <sequence>
        <any namespace="##other" processContents="lax"
            minOccurs="0" maxOccurs="unbounded"/>
    </sequence>
</complexType>

Although by definition the CIS coverage metadata item can contain literally “anything” and, hence, is outside the scope of the Standard, other standards utilize the metadata in a normative manner. One example is the EU INSPIRE WCS standard which is based on CIS coverages with an INSPIRE-specific metadata contents.

Tools are expected to add further information. For example, rasdaman can add information about the contributing footprints of the original input files to allow tracing back the origin of pixel data if needed.

The INSPIRE Elevation Grid metadata are embedded as shown below, following the above metadata schema while honoring the INSPIRE schema.

<gmlcov:metadata>
    <el-covmd:ElevationGridCoverageMetadata xsi:schemaLocation="http://inspire.ec.europa.eu/schemas/el-covmd/4.0 http://schema.datacove.eu/ElevationGridCoverageMetadata.xsd"
xmlns:el-covmd="http://inspire.ec.europa.eu/schemas/el-covmd/4.0" …>
         …
    </el-covmd:ElevationGridCoverageMetadata>
</gmlcov:metadata>

The following example extends this: three different, independent metadata records coexist in the coverage.

<GeneralGridCoverage xmlns="http://www.opengis.net/cis/1.1/gml" ...>
    <DomainSet> ...  </DomainSet>
    <RangeSet> ...  </RangeSet>
    <RangeType> ...  </RangeType>
    <Metadata>
        <el-covmd:ElevationGridCoverageMetadata
            xmlns:el-covmd=
            "http://inspire.ec.europa.eu/schemas/el-covmd/4.0" ...>
            ...
        </el-covmd:ElevationGridCoverageMetadata>
        <card4l:Card4lMetadata xmlns:card4l="..." ...>
            ...
        </card4l:Card4lMetadata>
        <special:MySpecialMetadata xmlns:special="..." ...>
            ...
        </special:MySpecialMetadata>
    </Metadata>
</GeneralGridCoverage>

How to combine various metadata items in this “bag” is not standardized, and hence tools will have varying issues in extracting specific information again. It would be helpful, at a minimum, to define a structure enabling the coexistence of metadata items without impact. This appears a relatively low-hanging fruit: by introducing a convention of parallel, independent sub-elements in the metadata slot, tools can provide and consume exactly those metadata the users are interested in, without getting confused by all the other slots that may be present. Simple path expressions would allow selecting such sections.

In rasdaman, for example, there are two such compartments already available. One is reserved for technical use by rasdaman, such as for metadata applying only to specific regions in a coverage. Another one is for metadata according to the EU INSPIRE standards.

A proposal is that OGC establishes a registry, maintained by the OGC Naming Authority (OGC-NA) and overseen by the Coverages SWG as is the case with other coverage-related registry information already, where extensions can be registered. Such an extension would consist of:

• a name not already registered;

• an XML namespace defining the structure; and

• any further pertinent information, such as description, examples, authoritative experts, etc.

Based on such registries, tools can extract the metadata of interest and be sure these are not “polluted” by other structures stored with the coverage, and as such ignore all other metadata items not of interest.

This has been discussed for XML, but it should be done in sync also for the JSON and RDF encodings which are also specified in the CIS standard.

B.6.2.1.  Interpolation

In Earth sciences, coverage data in general describe spatially extended, time-varying physical phenomena, corresponding to the mathematical notion of a “field” in physics [61]. Sensors observing natural phenomena do not collect data only at the exact point of a direct position. Rather, date are integrated over some neighborhood of the direct position and present the result as the value at this position. Generally, these observations are aggregated over some region and time. Therefore, it may make sense to invert this quantization to discrete points and derive values for in-between positions from the direct positions available in the coverage.

Through interpolation, range values can be obtained for coordinates within the domain of a coverage which are not direct positions, usually by combining the values of several direct positions in the neighborhood of the coordinate location under inspection. Technically, interpolation is applied during resampling in reprojection and scaling operations where new direct positions are constructed which need to get values assigned.

As a first requirement, the domain must be in a CRS which allows expressing “in between” coordinates. Cartesian coordinates, for example, do not fulfill this while geographic and temporal coordinates do.

While computing such an interpolated value, an additional requirement is that the range type must support derivation of values. Only nearest-neighbor interpolation is always possible as it retains one of the original values. Further common interpolations, such as linear, quadratic, or cubic require additionally a continuous range data type supporting the interpolation math.

Obviously, the justification for interpolation cannot simply be deduced from the data type. Just because some range type is an integer does not preclude that new values can be derived through integer arithmetics. And just because the domain is continuous does not mean that interpolation can be applied freely. One example of a very careful handling of interpolation is kriging[62]. Kriging is a method similar to regression analysis developed originally for geostatistics (in particular, mining where underground data tend to be particularly sparse), but also applied in many other domains like hydrogeology, remote sensing, and astronomy. See [63] for an in-depth presentation of spatial data interpolation.

B.6.2.2.  Image Pyramids

A special interoperability situation occurs in the context of image pyramids. Image pyramids are used by virtually every tool to speed up zoom access which is why pyramids are considered a form of processing but are actually an auxiliary data structure. When a user zooms out, then data get requested at a lower resolution to reduce client load and transfer times. For faster response the server caches lower-resolution versions of the data so that scaling loads and processes only from the next-closest higher resolution and not the high-volume base level. The problem is that while deriving the pyramid levels, interpolation may be applied. This raises the following concerns.

• With two scale steps being applied to the data, is the combined interpolation equivalent to a single scale operation, especially when considering numerical effects?

• During preparation of the pyramid when some interpolation gets applied, a choice must be made. If the user requests another interpolation method than the one used in pyramid building, then two different interpolation methods get applied to the data in turn, with unclear outcomes. Obviously, maintaining pyramid variants for different interpolations is unfeasible in face of the storage needs.

Obviously, this holds not only for a visual client’s zoom-out, but also if the server uses pyramids for scaling operations during any other processing. Generalizing this idea, a coverage should only allow interpolation methods during retrieval and processing which are compatible with the interpolation methods applied during generation of this coverage.

In the case of Landsat 8, chosen as a random example, it is found in [64] that the TIRS bands 10-11 are collected at 100 meters but get resampled to 30 meters to match the OLI multispectral bands. Actually, handling different bands in satellite images with different methods is quite common, given that several different instruments contribute bands. For Landsat 8, resampling is reported to be Cubic Convolution.

B.6.2.3.  Summarizability

It is well known in statistics that aggregation of values is meaningful only under particular conditions, and different types of aggregation are subject to different constraints. A large body of research is available from OLAP and Statistical Databases, including [65][66][67][68][69][70][76][77][78]. These constraints are given by the type of data (the coverage range type), but also by the particular axis (the coverage domain), and aggregation constraints typically vary across the axes. Finally, the type of aggregation (count, min, max, sum, avg as well as modes, percentiles, etc.) is relevant. Examples for such impact factors include the following.

• On classification data, maximum and minimum are meaningful, but neither sum (as well as difference) nor average is. Counting of occurrences is allowed.

• Addition of timestamps is not meaningful — adding years 2023+2023 would result in the year 4046 which is likely not an intended result. Rather, it is necessary to differentiate between timestamps (where addition is not meaningful) and time periods (which can be added up and added to timestamps).

• To obtain the car throughput in an image timeseries it is incorrect to add the number of cars recognized in each time slice because cars driving through the area of interest likely will appear in several such slices.

The term summarizability has been introduced by Lenz and Shoshani [65] to denote whether some particular type of aggregation is meaningful. Categorial data, for example, are not amenable to summing up or averaging, but there are more cases. The paper lists a series of examples, such as in socio-economic data, “number of accident deaths by month” is summarizable whereas “number of children attending school by month” is not summarizable.

Mazon et al. [66] provide a review of summarizability issues with a focus on Online Analytical Processing (OLAP). In [67] non-, semi-, and fully-additive measures are discussed, and beyond sums also averaging and rounding. Manifold approaches have been taken to ensure summarizability, including analytical [68], rule-based [69], and SQL integrity rules [70]. Inaccurate summary factors and the impact of inaccurate summary factors are treated in [76], providing practical rules on proper summation design. Altogether, there is a rich body of research and results available for business and statistical data.

Lenz and Shoshani [65] provide a categorization into flow (de-noting a cumulative effect over a period), stock (state at specific points in time), and value-per-unit, based on which they establish summation rules and compatibility matrices proving that these are necessary; no proof is provided for sufficiency. Kimball and Ross [77] adopt a slightly different perspective in distinguishing additive, semi-additive, and non-additive categories. On a side note, summarization does not preserve proportions [80], a fact which potentially can be relevant for model-based processing of EO data.

Niemi et al [78] combine these insights into the observation that a query can only produce meaningful results if (i) the aggregation operation is appropriate for the measure (in coverages: range values) and (ii) the measure is appropriate for the aggregation levels in the cube’s dimensions (in coverages: axes).

Niemi et al [78] further establish a new categorization of dimension types and multi-dimensional structures, derive a measure categorization, give formal definitions for summarizability types based on the relational model of data, and then construct provably correct rules for summarization.

In the context of EO, direct aggregation is concerned as well as indirect aggregation performed in “zonal” operations like reprojection and scaling, plus “global” operations like radiometric corrections.

In the field of Earth data such considerations are at a very preliminary stage. OGC Sensor Web Enablement (SWE) Common [79] distinguishes an eclectic mix of Count, CountRange, Category, CategoryRange, Quantity, QuantityRange, Time, TimeRange, Boolean, and Text which does not appear helpful in providing summarizability guidance to tools.

More research is required to ensure that results can be transferred from the work mentioned or to perform any necessary adjustment. For example, if the data comprises wind directions taken at time intervals and there is a need to determine the prevailing wind during some period of time, then it is not appropriate to compute the sum or the arithmetic mean, but it is appropriate to compute the mode [78].

A goal should be to define a framework for categorizing data with respect to the data’s summarizability which should allow finding provably correct constraints ideally derived automatically for some given coverage. In the SQL ecosystem, critical errors can be discovered automatically and, hence, avoided[81].

A first, admittedly simple and incomplete, step would be to reflect the well-known statistics categorization into categorical, ordinal, and numerical data by adding such an attribute into the coverage range type, plus corresponding protection of aggregation, interpolation, etc. With an increased understanding and better metadata such safeguarding might successively get increased.

B.6.2.4.  Dimension Hierarchies

On dimensions, aggregation often is driven by the conceptual model of the axis. For example, over time it is common to aggregate up to days, up to months, or up to years. Aggregations might consider all data of some chosen interval or regular parts thereof, such as “every February.” A roll-up from days to weeks in OLAP sales data is mathematically very similar to a map zoom-out by a factor 7, maybe with the particularity that OLAP data regularly are null at weekends. It is typical to have several successive abstraction levels, commonly called dimension hierarchies or concept hierarchies.

Looking more closely, the time axis reveals a particularity: As weeks and months do not align in the calendar, OLAP time hierarchies are not strict (i.e., tree-like), but directed acyclic graphs (DAGs) where users can decide to roll up over weeks or over months. Figure B.8 schematically shows two dimensional hierarchies.

Figure B.8 — Sample dimension hierarchies for geographic names and time

For Earth data, use of the temporal axis is likewise found. However, spatial x/y/z axes differ in the zoom semantics. As opposed to many other axis types, these are continuous in aggregation, or, expressed in common terms, allow zoom by any real-valued factor. (This is not to be confused with image pyramids where the predefined pyramid members only act as an internal acceleration and optimization mechanism, but do not serve any conceptual purpose and are not visible to users.)

Altogether, dimension hierarchies are powerful concepts for adding more semantics to multi-dimensional data (including summarizability constraints, as discussed before). It is certainly worth looking into solutions in other domains, such as ISO SQL OLAP [82] and Microsoft MDX [83]. The recommendation, therefore, is to add dimension hierarchies to coverage axes and exploit their semantics in the GeoDataCube queries of ISO 19123-3 [58].

B.6.2.5.  Validity and Reliability Masks

One way of characterizing validity of range values is to annotate every cell (pixel, voxel, etc.) with information about its validity which may consist of global metadata identical for all direct positions, or characterizations individual for direct positions.

For example, the SoilGrids service [84] offers an uncertainty layer which can be activated by the user for visual inspection (see Figure B.9).

Figure B.9 — SoilGrid uncertainty layer visualization

The l2gen tool [85][86] which processes Level 1 data to Level 2 generates on the side a quality bit vector per pixel with a series of relevant information, including the following.

• Pixel is over land/over sea

• Sunglint

• Reflectance exceeds threshold

• Observed radiance very high or saturated

• Sensor view zenith angle exceeds threshold

• Solar zenith exceeds threshold

• Probable straylight contamination

• Probable cloud or ice contamination

• Atmospheric correction failure

• Very low water-leaving radiance

• Navigation failure

• Navigation quality is suspects

In WCPS, such masking can be used to disregard “bad” pixels, say, in an aggregation, expressed as

for $s in ( MyScene ), $m in ( MySceneMask )
return count( $s * ($m.sunglint = 0) )

In case of a statistical mask with, say, percent values between 0 and 100 a threshold might get applied in a query:

for $s in ( MyScene ), $m in ( MyScenePercentage )
return max( $s * ($m > 70) )

As can be seen, there is a much wider range of possible irregularities than just clouds and cloud shadows, which again requires some remote sensing experience to decide about relevance of individual flags for a given analysis task. As a Holy Grail, the system could determine which quality filters have to be applied for a meaningful result (and whether the remaining values still are sufficient).

Investigation is recommended on adding quality measures (as described or different) to coverages and investigate on automatically evaluating such quality measures in WCPS, maybe using ontologies. However, Machine Learning (ML) is not recommended.

Figure B.10 — Sentinel-1 scenes delivered by ESA with different processing parameters applied

B.6.2.6.  Product Provisioning Coherence

Data provenance is important information for assessing the usefulness of a particular data set for a given purpose. In particular, the processing applied can have an impact as studies [87] convincingly show. Even changes in software versions or processing parameters can have critical consequences. An example is given in Figure B.10 where both times Sentinel-1 radar satellite data are shown: left with an acquisition (and, hence, ESA processing) time between June 2017 and May 2018, right with a Sentinel-1 patch from 2023 extracted from the ESA Copernicus archive. Neural networks trained on the left-hand side image generation by experience fail miserably when applied to data of the generation shown on the right-hand side.

In the FAIRiCUBE project, several data sets were harvested from the European Environmental Agency (EEA). While building up datacubes from the the file sets, several inconsistencies were noted. In some cases, resolution got enhanced at some time. However, the dataset was continued to be advertised as the same and with the same name. In land use classification data, at some point in time new classes were introduced (fortunately none were deleted). In one dataset, the data semantics changes over time: The imperviousness dataset for the year 2018 has percentage values from 0 — 100% whereas for 2006 it has class numbers indicating 0-20%, 20-40 %, etc. Sometimes the legend was observed to change over time, and it remains to be investigated deeper whether this hints at deeper semantics changes in the data.

From the usability perspective and the coverage homogeneity idea, such changes establish new products — tools running automatically over the dataset time extent can generate disastrous errors without a chance of recognizing and reporting this.

Changes are not captured by the coverage constituents currently, and generally are hard to capture. Basically, a provenance record would need to capture not only the tools used, but also the software version and the concrete parameters used. However, this still leaves open the interpretation of how the changes affect the product properties.

Generally, the recommendation is to define processing levels rigorously enough so that measures for variation introduced by the processing get provided alongside with the result, or ideally variations are ruled out in the first place. Towards this Holy Grail, Strobl suggests a revision of the ubiquitously used processing levels [93]. The suggestion distinguishes processing along two “dimensions” (not related to coverage dimensions), measurand and spatiotemporal, which together span a 2D matrix of combinations where each field defines a possible processing status. This allows sorting the usual processing levels into the matrix, relating them, and spotting combinations not addressed by processing levels hitherto.

From a coverage standards perspective, Strobl’s measurand dimension relates to the coverage range set whereas the spatiotemporal dimension relates to the coverage domain set. However, while the disentangling of ARD impact factors is a good way forward there is no clear requirements analysis which could justify the distinction made between combination-ready, fusion-ready, analysis-ready, and inference-ready.

B.6.2.7.  Numerical Effects

In the processing heavy EO world and the common use of floating-point numbers, a series of undesirable effects can occur which falsify results. Sometimes this influence is negligible, sometimes the result becomes completely wrong. It depends on the degree of deviation from the theoretically correct result, but also on the accuracy requirements of the problem to be solved.

One blatant example is defining null values as floating-point numbers. This is dangerous because (i) many numbers cannot be represented exactly (such as the beloved decimal-based fractions such as 0.1) and (ii) any operation applied can introduce inaccuracies. Stull [94] reports a simple example: add up 30 times 0.1 and you will get 2.9999993, rather than 3. This is with 32-bit floats; with 64-bit floats the result is 3.0000000000000013.

Now imagine a null value is defined as 3. Both computations would miss the null value. One may consider this example un-realistic, but it is felt that it is not: During import of data into a service it is common to preprocess from raw to say Level 2 (see discussion above), with all the numerical issues. Obviously, computations may accidentally distort the outcome and thereby generate a value that is treated as non-null, for example, in aggregations. Conversely, values which should count as non-null might accidentally get mapped to null.

This constitutes a huge and dangerous source of errors which often do not receive the attention required, among others, because data scientists rarely have an education in Numerical Mathematics as computer scientists usually get during undergraduate studies. In general, a rather naïve handling of floating-point numbers is often observed.

One of the insights Numerics teaches is interval arithmetics, used for floating-point numbers [95][96]. In this approach, numbers are considered inexact in the sense that the numbers do not have s single known value, but “smear” out over an interval of validity, in practice given by an interval. This can capture the behavior of floating-point numbers with finite accuracy. Error propagation rules can help estimating the ultimate accuracy of float computations [97].

While interval computations are significantly more involved (and require significant skills on the developer side) and drag down performance as compared to point Numerics, interval computations provide an automatic precision control which supports assessing processing chains, among others.

While tools generally do not implement this approach due to the above obstacles mentioned, these concepts might at least give hints on proper handling of finite-accuracy numbers. In rasdaman, for example, null intervals can be defined [98] which, when wisely used, can capture numerical deviations in the proximity of the conceptually foreseen null values. In the above example, a null interval around 3.0 might be defined as [2.9999993 : 3.0000007], based on the known processor accuracy. If a null value is defined as 3 then an equality comparison will fail in both cases.

The situation is further complicated through a parallel execution of code on heterogeneous hardware which can happen locally when CPU and GPU employ slightly divergent accuracies (and, hence, arithmetics), and even more so with distributed processing where, in the extreme case, supercomputer and edge devices jointly solve processing tasks.

Generally, for enhanced confidence, data tools should employ interval arithmetics so that the tools can provide error estimates. If provided by all tools along a processing chain, then the chain itself could be characterized in its overall error estimates. One starting point is shown with rasdaman where null intervals can be defined [98] and are strongly recommended on float and double range types.

B.6.3.1.  Service Quality Parameters

There are manifold general service quality parameters such as bandwidth, availability, and response time from coverage-specific quality parameters. Only the latter are of concern in this context which supports concentrating on the coverage structure by inspecting every normative component for accuracy and fitness considerations. In the following incomplete walk-through, qualitative and quantitative criteria for fitness of purpose of a given coverage for some applications are addressed. Further research should conduct a systematic review with ideally formal arguments for correctness and completeness.

• Domain

  • Match with the area of interest: are the requested dimensions present, such as a temporal axis for timeseries analysis?

  • Is spatiotemporal resolution sufficient for the purpose?

  • Are data provided in a CRS which supports transformation into one of the CRSs understood by the client?

• Range type

  • Are data in some unit that can be converted to a unit understandable by the client? For example, feet can be converted to meters, but flight (pressure) levels cannot easily be converted into hekto-Pascal (hPa).

  • Accuracy: Are the range values exact enough for the analysis to be performed so that the targeted results can be achieved? In the simplest case, this would become a parameter in the range type. However, accuracy may vary across the coverage extent, so other means might be appropriate such as accuracy masks. Note that this is related to null values (a null value has an accuracy of 0), but not identical (non-null values can vary wildly in their accuracy).

  • Plausibility: With what probability can the data be trusted? Again, this can be represented by a single value or, in the other extreme, by a twin coverage containing a probability for every range value. Measures can be binary or percents, among others (such as the floating-point accuracy or smaller).

• Range:

  • Are there sufficient non-null values in the area of interest, in all dimensions?

• Processing:

  • Are all aggregation and other functions applicable on the type of data (categorial, ordinal, numerical)?

  • What interpolation must be applied, if any? Is a suitable interpolation method admitted by the coverage?

The key question is always: is the given coverage good enough for answering the question the client has — where “good” can mean a variety of things, as discussed?

Which, however, prompts for a counterpart on the client side: What is the required minimum accuracy for a given task? This question has not been addressed adequately to the best of knowledge. More generally, this addresses service quality parameters which actually may have strong mutual dependencies. For example, JPEG image quality may be reduced to obtain smaller files which transmit faster. This quality reduction (in this case: by filtering out higher spectral frequencies, but also introducing artifacts sometimes) may or may not be acceptable.

Therefore, the development of APIs is suggested which, in parallel to formulating extraction and processing requests, enable clients to express quality criteria. Corresponding frameworks exist already in the field of cloud computing [99][100] including, for example, agent-based fuzzy approaches [101]. For semantic handling in the coverage context these need to be augmented with coverage-specific parameters, such as the ones discussed above. Optimal flexibility is given by providing a language rather than just a list of parameters so that conditions of any complexity can be built easily. For the evaluation of such quality requests against the quality offered by the service — essentially, a multi-attribute prioritization decision making problem — methods such as dual hesitant fuzzy graphs [102] may turn out promising.

B.6.3.2.  Coverage Fusion

As discussed in the Introduction of this Annex, data fusion is an excellent study field for ARD. Assume two georeferenced grid coverages, A and B. The goal is to obtain the mean square difference over all cells. In WCPS, this is expressed as:

for $a in (A), $b in (B)
return avg( pow( $a - $b, 2 ) )

If this is done as part of some large, unsupervised analysis process several facets must be investigated by the code (rather than a human), simply derived from the coverage definition. First, each coverage needs to undergo individual scrutiny as described in the previous subsection. Additionally, the following criteria are spotted; as before, rigorous research should establish a formally correct and complete criteria set.

• Domain

  • Do A and B share the same dimension and axes? Otherwise, the coordinate tuple is ambiguous.

  • Do both share the same CRS? If a reprojection is required, this could be performed automatically.

  • Do further constituents fit in a compatible way, such as: coverage type; categorial vs ordinal vs numerical; etc. Again, in some cases an automatic adjustment may be possible.

  • In case of resampling, are the interpolation methods of both coverages compatible?

• Range type

  • Through type casting, different numerical types (integer, float, complex) get adjusted automatically. This is part of the WCPS Standard already.

  • Determine proper treatment of null values (which may be different for A and B).

  • Verify the units of measure for compatibility between A and B. If there is an opportunity for automatic unit conversion, apply that. Further investigation should determine how much of this inter-coverage harmonization can be automated.

B.6.3.3.  ML

From the perspective of this discussion, model-based prediction through Neural Networks constitutes just another “processing” request to a service. For example, in the UDF-enhanced WCPS of rasdaman a prediction can be invoked as a query function which internally the server, as part of the query processing sends over to pytorch. The following example takes some user-specified patch from a datacube built from Sentinel-2 optical satellite data, and applies a crop model to it.

for $c in (Sentinel_2a), $m in (CropModel)
return encode( nn.predict( $c[…], $m ),“tiff“ )

Technically, function predict() from package nn (Neural Networks) receives the area of interest together with the model to be applied and passes both to pytorch. The result returned is processed further by rasdaman and in this case encoded as TIFF.

As it turns out in experiments as part of the AI-Cube [108] and FAIRiCUBE [109] projects, such pretrained models are extremely susceptible to even small deviations from the training data in practically all parameters: different resolution, different radiometry, different landscape, etc. Therefore, a current research topic is to investigate how to “fence” models in a way that the server can predict the result quality and possibly warn users or even refuse to apply a model in a given situation, such as “applying a fire forecast over the ocean does not make sense.” Initial-stage considerations on criteria for such an assessment include the following.

• Landscape type: As a simple example, a seashore and water areas boundary polygon might automatically be applied on areas of interest when models classified as land-based are invoked.

• Temporal applicability: For example, a crop damage estimator should be applied at a time of year when normally crop grows in the area of interest.

• Data-inherent properties: This might address freeness of clouds and cloud shadows, data histogram, or other suitable statistics.

Obviously, to this end models also need to carry quality information, specifically: a description of applicability. This addresses on the one hand the input data for the model analysis, but also the quality parameters imposed by the client. It is understood that such research is in its infancy, although highly important when model-based prediction occurs somewhere in the middle of automated decision-making pipelines, without human plausibility checking.