B.3.1.  Syntax

aligned(8) class TiledImageConfigurationBox
extends ItemFullProperty('tilC', version=0, flags) {

  unsigned int(32) tile_width;
  unsigned int(32) tile_height;

  unsigned int(32) tile_compression_type;

  unsigned int(8) number_of_extra_dimensions;
  for (int i=0; i<number_of_extra_dimensions; i++) {
    unsigned int(32) dimension_size[i];
  }

  unsigned int(8) number_of_tile_properties;
  for (int i=0; i<number_of_tile_properties; i++) {
    ItemProperty tile_image_property[i];
  }
}

Listing B.1 — Syntax for TiledImageConfigurationBox

B.3.2.  Semantics

tile_width, tile_height is the size of a single tile. All tiles have the same size. Tiles at the right or bottom border may extend beyond the total image size. This is consistent with the definition of a tessellation type “square” in the OGC Abstract Specification Topic 22 — Core Tiling Conceptual and Logical Models for 2D Euclidean Space, OGC 19-014r3.

tile_item_type specifies the image item type used for all the individual tile images. tile_item_type is one of the possible four-character types of ordinary image items (e.g. hvc1 for h265 compression or j2k1 for JPEG2000).

number_of_extra_dimensions specifies the number of dimensions of the N-dimensional image as number_of_extra_dimensions = N - 2. A 2D image has number_of_extra_dimensions=0.

dimension_size[i] specifies the size of dimension i+2 of the N-dimensional image. The size of the first two dimensions is obtained from the mandatory ispe item property.

number_of_tile_properties specifies the number of tile properties stored in the tilC.

tile_image_property[] are the image item properties used when decoding a tile image. This includes at least all mandatory item properties for an image item of type tile_item_type with the exception of ispe. If tile_image_property[] does not contain an ispe property box, the decoder shall synthesize an ispe property with tile_width and tile_height as the size.

OffsetFieldLength = OFFS_LEN[flags & 0x03] defines the number of bits used to store the offset to the image data of a specific tile. OFFS_LEN[] = [ 32, 40, 48, 64 ]

SizeFieldLength = SIZE_LEN[(flags>>2) & 0x03] defines the number of bits used to store the length of the image data of a specific tile. SIZE_LEN[] = [ 0, 24, 32, 64 ]

(flags & 0x10) is a hint to a decoder whether the compressed tile image data is stored consecutively in sequential order.

The four different offset pointer sizes defined the maximum tili image file sizes as shown in Table B.1

Table B.1 — Maximum image file size depending on pointer lengths.

The four different tile size field lengths define the maximum compressed tile sizes as shown in Table B.2

Table B.2 — Maximum compressed tile size depending on size field length.

B.3.3.  tili Item Data

The item data consists of an offset pointer table TiledImageOffsetTable, followed by the compressed image data.

The number of tile offsets stored in the table (NumTiles) is computed by:

  TileColumns = (ispe_width + tile_width -1)/tile_width;
  TileRows    = (ispe_height + tile_height -1)/tile_height;

  NumTiles = TileColumns * TileRows;
  for (i=0; i<number_of_extra_dimensions; i++) {
    NumTiles = NumTiles * dimension_size[i];
  }

Listing B.2

ispe_width and ispe_height is the total image size as specified in the mandatory ispe item property.

aligned(8) class TiledImageOffsetTable {

  for (int i=0; i < NumTiles ; i++) {
    unsigned int(OffsetFieldLength) tile_start_offset[i];
    unsigned int(SizeFieldLength) tile_size[i];         // note: not present if SizeFieldLength==0
  }
}
// ... followed by compressed tile data ...

Listing B.3

B.3.4.  Semantics

tile_start_offset[i] points to the start of the compressed data of the tile. The position is given relative to the start of the TiledImageOffsetTable data. Note that this is not a file offset, but an offset into the item’s data that can potentially span several iloc extents. If a tile is not coded and the corresponding image area is undefined, the tile_start_offset[i] shall be 0. If a tile is not coded, but the displayed image should be taken from a lower-resolution layer (in a pymd stack), tile_start_offset[i] shall be 1.

tile_size[i] (if present) indicates the number of bytes of the coded tile bitstream.

The entries in the offset table are ordered in row-major sequence. In other words, for a 2D image, they are indexed as [y][x], a three dimensional volumetric image (extra dimension z) would be indexed as [z][y][x].

The compressed tiles data may be stored in the file in arbitrary order, i.e. the tile_start_offset[] are not necessarily in increasing order.

If the tile_size[i] variables are not present, the decoder has to infer them from the tile_start_offset[] values. For the case that the tiles are stored in sequential order ( flags & 0x10 == 0x10), the tile_size[i] can be computed as tile_start_offset[i+1] - tile_start_offset[i] with the exception of the last tile, which extends until the end of the data. If the tiles are not stored in sequential order, the decoder first must sort the tile start offsets before it can again compute the size from the difference to the next tile start. Note that in this case, the decoder cannot read the offset table on-demand. Thus, storing the tile sizes is advised in this case if on-demand access of the tile offsets is desired. Multiple tiles can use the same tile_start_offset to reference a similar image content. This case must be taken care of when computing the tile sizes.

The file structure of a tili image item is illustrated in Figure B.1.

Figure B.1 — File structure of a tili image item. Picture taken from m70021:Tili proposal, drawn by Joe Stufflebeam.

B.4.1.  Multi-Dimensional Data

With the extra_dimensions variable in the tilC configuration header, it is possible to store not only 2D image data in a tili image item but also organize the tiles into datasets with additional dimensions.

For example, one extra dimension could be used to stack the 2D tiled images into a 3D voxel image volume.

Another example would be to use the extra dimension to address different frequency bands in a multispectral image. This is illustrated in Figure B.2.

Figure B.2 — File structure of a tili image item. Picture taken from m70021:Tili proposal, drawn by Joe Stufflebeam.

This scheme can be extended to even more extra dimensions. As an example, when two extra dimensions i and k are used, the tiles in the TiledImageOffsetTable are indexed by [k][i][y][x].

D.1.1.  Two approaches for modular development

GDAL and Apache SIS have different development strategies. GDAL is a translation layer between a unique API and various independent libraries. This approach allows GDAL to support a large number of formats and allows a large part of the work to be delegated to other teams such as the libheif developers. By contrast, Apache SIS supports a smaller number of formats but often implements these formats itself rather than delegating to an external library. This approach looks like a much larger effort, duplication of work and risk of bugs. But actually, GeoTIFF, GeoHEIF and many other chunked formats have many fundamental principles in common. The formulas for translating areas of interest and subsampling into tile indices, scanline strides, pixel strides and band indices are the same, even if the way to use this information is different. Decompression algorithms such as “deflate” and JPEG operate on arbitrary streams of bytes, and therefore are either already format-independent, or can be made format-independent in the cases of algorithms originally developed specifically for GeoTIFF. For example, the algorithms for caching tiles and assembling them in larger images are also the same.

D.1.2.  Reusing a code baseline for common functionalities

By building a GeoHEIF reader on top of a shared code base, the development of a prototype supporting tiling was relatively quick. The majority of the development effort is in parsing metadata, which, in GeoHEIF case, are numerous.

E.1.1.  Why “pixel is point/area” are not synonymous of “pixel center/corner”

Consider the process of resampling a source raster in one CRS into a destination raster in another CRS. Let T be the coordinate operation for the whole chain of transforms, from pixel coordinates in the source raster to pixel coordinates they should have in the destination raster, with a map projection in the middle if desired. The resampling process will iterate over all pixels to fill in the destination raster and, for each destination pixel, use the inverse of the coordinate operation (T⁻¹) for computing the pixel coordinates where to look for a value in the source raster. The destination pixel coordinates will typically be integers (represented in the figure below by circles in pixel centers), since the process iterates over these pixels. But the corresponding source pixel coordinates (represented in the figure below by arrows) can be anywhere. They will usually be non-integers, thus requiring interpolations between two or more source pixel values.

Figure E.3 — Resampling a raster with a coordinate operation T from source pixels to destination pixels

Because bilinear or bicubic interpolations are computed between points, the resampling process needs to choose which point in cell will be used. Should T_(i,j) transforms from pixel center to pixel center, or from pixel corner to pixel corner? For the “pixel is point” footprint, mapping pixel centers seems a natural choice. However, also for the “pixel is area” footprint, mapping pixel centers is desired. If a measurement is the average value of the whole pixel area, then when the area of a destination pixel covers the areas of two source pixels, the destination pixel should be the average of the two source pixels. This is illustrated on the left side of the figure below, with the yellow color representing the average of the green and red colors.

Figure E.4 — Linear interpolation of pixel values with different "pixel center/corner" conventions

Computing the average of two values is mathematically equivalent to a linear interpolation evaluated midway between these two values. Therefore, mapping pixel centers produces the desired effect with linear interpolation for both “pixel is point” and “pixel is area” measurement footprints. If the transform was mapping pixel corners instead, as illustrated on the right side of above figure, the interpolation result would be the value of the leftmost pixel instead of the average value.

A resampling process may also need to estimate the real world bounding box of the result. For the sake of simplicity, consider pixels with a resolution of 1° in longitude. If the image width is 4 pixels, the real world width may be 4° or 3° of longitude depending on which side of Figure E.1 is applicable. Lets assume a coordinate operation scaling the resolution to 2° of longitude, as illustrated in Figure E.4. When using the “pixel corner” convention (i.e., the strategy illustrated on the left side of Figure E.2), the green box on the left side of the figure below is transformed to the green box on the right side. The result is 2 pixels in width and is aligned with the new grid. But when using the “pixel center” convention, the yellow box on the left side become the yellow box on the right side. The result has a width of 1.5 pixels and is not aligned with the grid, i.e. the bounding box starts at -0.25 pixel relative to the new pixel centers (red circles).

Figure E.5 — Bounding box transformation using "pixel corner" (green) and "pixel center" (yellow) conventions

It may be argued that for a “pixel is point” footprint, the yellow bounding boxes are physically more correct. However, for some purposes such as building an image pyramid, the green bounding boxes have more convenient properties. Therefore, the use of a “pixel corner” convention can be legitimate even with a “pixel is point” footprint.

E.1.2.  Ambiguities in axis order and directions

In the OGC API — Common (draft) and CIS models, the linear relationship from pixel coordinates (i, j, k, m) to real world coordinates (x, y, z, t) can be as below (lower_grid is zero in the OGC API case):

    x = lower_world[0] + (i − lower_grid[0]) × resolution[0]
    y = upper_world[1] − (j − lower_grid[1]) × resolution[1] (assuming inclusive upper bound)
    z = lower_world[2] + (k − lower_grid[2]) × resolution[2]
    t = lower_world[3] + (m − lower_grid[3]) × resolution[3]

Typically, only the two first dimensions are used in above list. This approach is generalizable to any number of dimensions, ignoring the ambiguity about whether other special cases like y may exist. However, this approach assumes that the real world coordinates are in the same order as the grid coordinates. This assumption is problematic and has been a great source of confusion for decades. This is because screen coordinates and geographic coordinates have inherited different conventions from their history. Most raster producers and users want the grid (image) axes organized in such a way that the image appears with a familiar orientation (e.g., north toward up) when shown in a web browser or other non-geospatial software. Since screen coordinate systems traditionally use (right, bottom) axis orientations in that order, the corresponding CRS axis orientations would be (east, south). However, many geographic CRS use (north, east) axis orientations instead (there is also some west-orientated and south-orientated geographic CRS). Not only the axis order differs, but the north-south orientation also differs. The “georeferencing by bounding boxes” approach usually tries to address this mismatch with hard-coded rules.

E.2.1.  Mapping grid axes to CRS axes

This section provides a step-by-step guide for defining an affine conversion in replacement of the “georeferencing by bounding boxes” approach. For the sake of simplicity, the example used in this section is a two-dimensional raster written in the GeoHEIF format with (latitude, longitude) CRS axes and with the (0,0) pixel coordinates located in the upper-left corner. However, the steps described here are generalizable to any number of dimensions, any grid corner location, any axis order and orientation, and any raster formats capable to store affine transforms.

The simple case shown in this section goes as follow:

The affine conversion can be built in a two-steps process as follow.

AffineTransform gridToCRS = ...;    // Using "pixel corner" convention.
gridToCRS.translate(0.5, 0.5);      // Now using "pixel center" convention.

E.2.1.4.  GML in JPEG2000 case

GML in JPEG2000 (derived from ISO 19123-2) splits the affine transform into two components: the offset vectors and the origin. These two components can easily be assembled into the matrix form used by this annex:

 

Figure E.6 — Offset vectors and origin in an affine transform

That “offset vectors + origin” way to specify the conversion is as good as the matrix form used in this annex. The separation between offset vectors and origin is sometime convenient in explanatory texts. For example, for saying that there is a strict correspondence between grid axes and CRS axes (no shear, rotation or axis swapping, but potentially axis direction reversal), the two following sentences are equivalent:

• The matrix is diagonal, ignoring the last column.

• The offset vectors part of the matrix is diagonal.

For specifying that axis swapping are authorized but not shears or rotations (except by multiple of 90°), the two following sentences are equivalent:

• Each matrix row and each matrix column shall contain exactly one non-zero value, ignoring the last column.

• Each offset vector shall contain exactly one non-zero value, and the index of that non-zero value shall be different for all offset vectors.

With the above rule, the interpretation of the offset vectors part of the matrix is as follows. For any row j having exactly one non-zero value in column i, that value is the scale factor no matter in which column it appears. Index j is the CRS axis index and i is the grid axis index. The scale factor is multiplied by the grid coordinate at index i and the result is stored in the CRS coordinate at index j. When the scale factor is on the diagonal, the grid and CRS axes are at the same position (j = i), i.e. these axes are aligned.

int dimension = 2;          // Could be anything, this algorithm is generic.
double[][] matrixInGeoHEIF = ...;
double[][] matrixInGeoTIFF = new double[dimension][];
for (int inEffectiveCRS = 0; inEffectiveCRS < dimension; inEffectiveCRS++) {
    int inRealCRS = ...;    // Find the index in the real CRS having the same axis direction.
    matrixInGeoTIFF[inEffectiveCRS] = matrixInGeoHEIF[inRealCRS];
}

Data axis to CRS axis mapping: 1,2
GeoTransform =
  90, 0, -0.3515625
  -180, 0.3515625, 0

GeoHEIF::InitializeCRS() {
    OGRSpatialReference crs = ...;          // Parse the CRS from the GeoHEIF file.
    crs.setAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
}

GeoHEIF::GetGeoTransform(double *padfTransform) {
    OGRSpatialReference crs  = ...;         // Above CRS with OAMS_TRADITIONAL_GIS_ORDER.
    int[]      dataAxisToSRS = crs.GetDataAxisToSRSAxisMapping()
    int[]      rasterSize    = ...;         // Raster width and height in pixels.
    double[][] matrixInHEIF  = ...;         // The matrix read from GeoHEIF.
    double[][] matrixInGDAL  = new double[dimension][];
    for (int j=0; j<dimension; j++) {
        int dataIndex = abs(dataAxisToSRS[j]) - 1;
        matrixInGDAL[dataIndex] = matrixInHEIF[j];      // Take a matrix row.
        /*
         * By GDAL convention, a negative value means that the axis direction is reversed.
         * This implies that the scale factor shall be negated and the translation changed
         * from minimum CRS value to maximum CRS value (if the scale factor was positive).
         */
        if (dataAxisToSRS[j] < 0) {
            double span = 0;
            for (int i=0; i<dimension; i++) {
                double coefficient = matrixInGDAL[dataIndex][i];
                matrixInGDAL[dataIndex][i] = -coefficient;
                span += coefficient * rasterSize[i];
            }
            matrixInGDAL[dataIndex][dimension] += span;
        }
    }
    // Copy values to `padfTransform`, keeping in mind that GDAL puts translation first.
}

E.2.2.  Transforms chain from raster to screen

The affine conversion built in Annex E.2.1 is describing the conversion from raster coordinates to CRS coordinates, not to screen coordinates. In this section, the above-cited conversion is called the Grid to CRS conversion. For visualization, another step from CRS coordinates to screen coordinates is required. That second step can also be expressed as an affine transform and is called the CRS to Display conversion in this section. Those two conversions have clearly separated tasks:

Note that the CRS to Display conversion is the same for all data to be rendered on screen: The other rasters, the geometries, and all feature as long as they are using the same CRS. Note also that this conversion is always needed no matter how the data are georectified (by bounding boxes or by affine transforms). This is not a complication introduced by the “georeferencing by affine transforms” approach. The only difference is that this section will describe that conversion more rigorously using affine transforms.

For the first part of this section, the “CRS” term in “Grid to CRS” and in “CRS to Display” must refer to the exact same CRS, axis order included. If that CRS uses (latitude, longitude) axis order, then the “CRS to Display” conversion also needs to swap axes to the (x, y) order. The case where the CRS are not the same will be discussed in Annex E.2.2.3.

RenderedImage   raster       = ...;     // Read from GeoHEIF file.
AffineTransform gridToCRS    = ...;     // Read from GeoHEIF metadata.
AffineTransform crsToDisplay = ...;     // Computed by the viewer.
Graphics2D      output       = ...;     // Typically provided by widget.

output.transform(crsToDisplay);
output.drawRenderedImage(gridToCRS, raster);    // Java2D does matrix multiplication for us.

// Example of a geometry in geographic coordinates (EPSG:4326).
var shape = new GeneralPath();
shape.moveTo( 30, 50);      // 30°N, 50°E
shape.lineTo(-20, 60);      // 20°S, 60°E
shape.lineTo(-10, -5);      // 10°S,  5°W
output.draw(shape);         // The "CRS to Display" transform applies also to this shape.

E.3.1.  Transition between map projections

When reprojecting the raster data to another CRS, Apache SIS updates the “CRS to display” transform for the change of CRS in such a way that the visual change is imperceptible at the location of the mouse cursor. To demonstrate:

Note that the feature where the right-click occurred seems unchanged. If the raster data covers a small area, it may look like nothing happened at all. On raster datasets large enough, the change become more and more apparent as the distance from the right-click position increases.

While the current version of the JavaFX application does not yet show the following statement, Apache SIS is capable of following the displacement of a vehicle with automatic changes of UTM zone. This could also show automatic changes from UTM to stereographic projection when the vehicle enters polar regions, and keep the changes visually unnoticeable by users — at least in the vehicle’s vicinity. This feature requires the support of Jacobian matrices (not discussed in this annex), which can be seen as a complement to affine transforms for non-linear cases. The above scenario works for any map projection in which Jacobian matrix support has been implemented. Jacobian matrices were introduced (under a different name) more than twenty years ago in OGC 01-009.

E.3.2.  Positional accuracy shown in coordinate values

The bottom right corner of the image below shows the geographic or projected coordinates of the mouse cursor’s position. Each coordinate value uses, independently of other values in the coordinate tuple, the minimal number of digits needed for distinguishing two screen pixels at the cursor’s position. This number of digits varies not only with the zoom level, but also simply by moving the mouse to different regions of the displayed raster dataset. To demonstrate:

This is done by inspecting the scale factors (a little bit more than only the diagonal) of the “CRS to display” matrix. If the chain of operations involves a datum change, Apache SIS also compares the scale factors with the operation positional accuracy and appends a text such as “± 100 metres” when the pixel size become smaller than the operation accuracy.

Figure E.7 — Demo application showing positional accuracy

E.3.3.  Synchronized user’s action between windows

Another affine transform usage demonstrated by the Apache SIS viewer is the synchronization of the geographic areas shown in different windows. When many windows are showing the same geographic area of different rasters, and when the user changes the viewed area in one window, it is easy to propagate that change to all windows if the same “real world” bounding box and the same orientation (rotation) is applied everywhere. However, when the windows are showing data in different geographic areas (for example, one window showing an area at the south of another window), potentially using different orientations or even map projections, synchronizing the windows become more challenging. It is possible to perform those more advanced synchronizations by expressing every user’s action, such as zoom or pan, in one window as an affine transform. That affine transform can be converted from pixel units of one window to real world units, then to pixel units of another window. The result can be applied to the other window even if the two windows are showing different data at different resolutions with different zoom levels, geographic areas, rotations or map projections. To demonstrate:

This synchronization is done by expressing all user’s actions as affine transforms, converting these transforms with matrix operations (inversions and multiplications), then concatenating the converted actions to the affine transform of another window.

Figure E.8 — User’s action propagated to different windows

F.2.1.  “NaNs are toxic”

NaNs change everything they touch to NaN. But this is exactly the desired feature. For example, in (x + y) / 2, if y is unknown, then the result is unknown. Without NaNs, the classical approach is to use Sentinel’s “no data” values such as -9999 like below:

double average(double x, double y) {
    if (x == -9999 || y == -9999) {
        return -9999;
    } else {
        return (x + y) / 2;
    }
}

Listing F.2 — Example of code handling a Sentinel’s "no data" value

NaNs make the above verifications unnecessary. This is not only a convenience, but also (most important) a safety. In complex applications, it is very likely that some checks will be forgotten, or that the developer will wrongly assume that -9999 will never reach some point. This risk can be compared with buffer overflows caused by missing bound checks, which are still a major cause of bugs and security breaches despite decades of developers effort. Those problems can go unnoticed for some time before they cause harm.

Example 1: When computing the average value of {23, 33, -9999, 30, 34} where -9999 is a Sentinel’s “no data” value that the developer forgot to handle, the result is -1975.8. At this point, the developer may suspect that something is wrong, but the software will probably not. However, if the average is computed over 1000 random values in the [10 … 50] range, then the corrupted result is close to 20 (by comparison, the normal value is close to 30). This corrupted result is well inside the range of valid values, potentially making the error unnoticeable even for the developer.

Example 2 (could be a demonstration): The demo programs check for missing values in raster data. However, it does not check for missing values in the coordinates of the points where to interpolate. If those coordinates were latitudes and longitudes, and if the calculation involved a map projection where all latitude values were given to sin, cos or tan trigonometric functions (this is a common scenario in map projections), then -9999° would be fully equivalent to 81°. It can make corrupted values totally indistinguishable from valid values.

Such floating-point errors should not be taken lightly since, for example, the first Ariane V rocket exploded because of a floating-point overflow. The use of NaN values reduces at least the risk that corruption caused by unchecked “no data” values stay unnoticed. Even if a library performs very carefully all “no data” checks, if the library and the caller using the library do not agree on the same Sentinel’s “no data” values (for example because of a misconfiguration), the result is the same as if the library performed no check at all.

F.2.2.  “There is only one NaN, while we need many no-data values”

This idea is false. The IEEE 754 float type has 8,388,608 distinct quiet NaN values (ignoring signaling NaNs, see below for notes on the difference) while the double type has over 10¹⁵ distinct quiet NaNs. When the bit pattern of a float is interpreted as an int, all values in the range 0x7fc00000 to 0x7fffffff inclusive are NaN values. This is not the only range of NaN values, but this is the most convenient range, and the one used in the OGC Testbed 20 demonstration. See Annex F.3.1 for an example of definition of NaN values in that range.

Verification: for demonstrating the availability of many distinct NaN values, the OGC Testbed 20 demonstration ( Annex F.5) used four of values:

To demonstrate that NaN values can be analyzed, the bilinear interpolation arbitrarily applies the following rule: If exactly one value needed for the interpolation is missing, the interpolation result is missing for the same reason. If two or more needed values are missing, the interpolation result is missing for the reason having the highest priority in the following order (highest priority first): NO_PASS, LAND, CLOUD and UNKNOWN. Executions and comparisons against the precomputed expected values show that the different NaN values are correctly distinguished in all tested languages.

F.2.3.  “Not the right semantic (we want null, unknown, …)”

NaN means “Not A Number”. There is no semantic associated to why a value is NaN. It can be any reason of the user’s choice. With more than 4 millions distinct NaN values available, there is plenty of room for assigning whatever semantic the user wants to different NaN values. Assigning a “no data” meaning to, for example, -9999, which is a number contrary to NaN, is not semantically better.

NaN cannot represent directly a JavaScript null or unknown value (neither can real values such as -9999). However, null and unknown semantics are a layer on top of the IEEE 754 values. For example, in the Java and C/C++ languages, null values cannot be assigned to primitive types. In these languages, null values can only be assigned to wrappers (Java) or pointers (C/C++). For example, a java.lang.Float (in Java) or a float* (in C/C++) can be null, but a float cannot. In JavaScript, the fact that null and unknown apply to all kinds of objects is an indication that they are handled at a level other than a primitive type.

F.2.4.  “The payloads that distinguish NaN values may be lost”

This argument is based on the fact that IEEE 754 recommends, but does not mandate, NaN’s payload propagation. This issue is discussed in Annex F.2.5, but does not apply to loading or saving values from/to files. When reading an array of float values, the operating system is reading bytes. If data are compressed with a lossless algorithm, decompression algorithms such as GZIP usually operate on bytes. Next, the conversion from big-endian byte order to little-endian or conversely is applied on bytes. Only after all these steps have been completed is the final result interpreted as an array of floating-point values. Casting a byte* pointer to a float* pointer in C/C++ does not rewrite the array, it only changes the way that the computer interprets the bit patterns that are stored in the array. Even after the cast, if a float array needs to be copied, the memcpy (C/C++) or System.arraycopy (Java) functions copy bytes without doing any interpretation or conversion of the bit patterns. Nowhere does the Floating-Point Unit (FPU) need to be involved. Therefore, there is no way that the NaN’s payload can be lost with the above steps.

Verification: The demonstration ( Annex F.5) provides the same raster data in two versions:

The demonstration reads the binary files, changes the byte order if needed, then interprets the result as floating-point values. See for example the loadRaster() method in the TestCase.cpp file, which ends with return reinterpret_cast<float*>(bytes). In the Java case, the read operation involves a copy from the byte[] array to a newly allocated float[] array, but execution shows that no payload is lost.

F.2.5.  “Payload propagation is implementation-dependent”

When a quiet NaN value is used in an arithmetic operation, IEEE 754 recommends that the result is the same NaN as the operand. However, this is not mandatory and the actual behavior may vary depending on the platform, language, or applied optimizations. Furthermore, the situation is even more uncertain when two or more operands are different NaN values. Which one takes precedence? Empirical tests suggest that, at least with Java 22 on Intel® Core™ i7-8750H, the leftmost operand takes precedence in additions and multiplications, while the rightmost operand takes precedence in subtractions and divisions. However, despite this uncertainty, users still have the guarantee that the result will be some NaN value. The fact that users do not know for sure which NaN values they got is not worse than the approach using Sentinel’s “no data” values. For example, when computing the average value of {2, 8, 9999, 3} where 9999 represents a missing value, if developers forget to check for missing values, they will get 2503. Subsequent operations consuming that value are likely to fail to recognize 2503 as a missing value, leading to unpredictable consequences. By contrast, when using NaNs, developers get some NaN result even if they forgot to check for missing values. Even if the users are uncertain about which NaN value was obtained, this situation is still better than an accidental value indistinguishable from real values. The preceding example only computed the average of 4 values, but the more values to average there are, the more difficult it becomes to notice Sentinel’s “no data” values polluting the computation.

Verification: The demonstration ( Annex F.5) using Sentinel’s “no data” values needs to check for missing values before using them in calculations. The same is true (in a more lenient way) for applications using NaN values if they wish to distinguish between the different NaN payloads. However, users of NaN values have two additional alternatives that users of Sentinel’s “no data” values do not have:

The TestNodata demonstration had to check for Sentinel values before any computation. This was no choice. The TestNaN demonstration could have done the same check but instead arbitrarily chooses to check for NaN payloads after calculations. This is an arbitrary choice, TestNaN could have been kept in a form almost identical to TestNodata. This is done for illustrating the flexibility mentioned in previously made points. This strategy can provide a performance gain when the majority of the values are not missing.

F.3.1.  “No data” values declaration

The TestNodata class uses the following Sentinel values for missing data:

/**
 * Sentinel value for missing data. A value may be missing for different reasons,
 * which are identified by different Sentinel values.
 */
static final float UNKNOWN = 10000,
                   CLOUD   = 10001,
                   LAND    = 10002,
                   NO_PASS = 10003;

/**
 * The threshold used for deciding if a value should be considered as a missing value.
 * Any value greater than this threshold will be considered a missing value.
 */
static final float MISSING_VALUE_THRESHOLD = UNKNOWN;

Listing F.3 — Example of distinct Sentinel’s "no data" value definitions

Note that the MISSING_VALUE_THRESHOLD value is arbitrary and data-dependent, as it must be a value not used by real values. Changing the threshold can require changes in the code, for example if the threshold becomes a value smaller than all the valid values instead of greater than them. By contrast, an approach based on NaN values does not need such an arbitrary choice. The TestNaN class uses the following NaN values for missing data:

/**
 * Value of the first positive quiet NaN.
 */
private static final int FIRST_QUIET_NAN = 0x7FC00000;

/**
 * NaN bit pattern for missing data. A value may be missing for different reasons,
 * which are identified by each of the different NaN values.
 */
static final int UNKNOWN = FIRST_QUIET_NAN,      // This is the default NaN value in Java.
                 CLOUD   = FIRST_QUIET_NAN + 1,
                 LAND    = FIRST_QUIET_NAN + 2,
                 NO_PASS = FIRST_QUIET_NAN + 3;

Listing F.4 — Example of distinct NaN value definitions

Alternatively, C/C++ developers can use the std::nanf(char*) function instead. This function is defined by the C++ 11 standard and an usage example is given in Annex F.4.

F.3.2.  “No data” values checks

The TestNodata class uses the following code for checking for missing values. In this case, this check must be done before using the values in formulas. The use of max is an optimization based on the fact that, in this test, the reasons why a value is considered missing are sorted in order of precedence (see Annex F.5.2 for more details). Production codes may be more complex if they cannot rely on this assumption.

float v00, v01, v10, v11 = ...;   // The raster data to interpolate.
float missingValueReason = max(
        max(v00, v01),
        max(v10, v11));

double result;
if (missingValueReason >= MISSING_VALUE_THRESHOLD) {
    // At least one value is missing.
    result = missingValueReason;
} else {
    // All values are valid. Interpolate now.
    result = ...;
}

Listing F.5 — Example of checking Sentinel’s "no data" values before interpolation

The TestNaN class uses the following code for checking for missing values. In this case, contrary to TestNodata, developers are free to check before or after any interpolation is computed. This example arbitrarily performs the check after the interpolation. The use of max below is the same trick as the one used above. Production codes may be more complex for the same reasons but have no reason to be more complex than the code using “no data” Sentinel values.

float v00, v01, v10, v11 = ...;   // The raster data to interpolate.
double result = ...;              // Interpolate regardless if values are valid.
if (isNaN(result)) {
    // At least one value is missing.
    // This block could be completely omitted if the developer
    // does not care about the reason why a value is missing.
    int missingValueReason = max(
            max(floatToRawIntBits(v00), floatToRawIntBits(v01)),
            max(floatToRawIntBits(v10), floatToRawIntBits(v11)));

    // Result was already NaN, but replace with another NaN for the selected reason.
    result = intBitsToFloat(missingValueReason);
}

Listing F.6 — Example of checking NaN values after interpolation

F.3.3.  Multiple missing value reasons for the same pixel

Handling NaN values as bit patterns allows a more advanced use case which is not possible (or at least not as efficiently) with “no data” Sentinel values. Instead of associating the missing value reasons to different NaN values (over 2 million possibilities), users can associate these reasons to different bits. The NaN payload has 21 bits (ignoring the signaling NaN bit and the sign bit), which gives room for 21 different reasons (or 22 if the sign bit is also used). The TestNaN code declaring NaN values can be replaced by the following:

static final int UNKNOWN = FIRST_QUIET_NAN,     // No bit set. This is the default NaN value in Java.
                 CLOUD   = FIRST_QUIET_NAN | 1,
                 LAND    = FIRST_QUIET_NAN | 2,
                 NO_PASS = FIRST_QUIET_NAN | 4,
                 OTHER   = FIRST_QUIET_NAN | 8;

Listing F.7 — Example of NaN definitions allowing multiple reasons to be packed in a single value

And the code computing the missing value reasons (previously using max) become the following:

if (isNaN(result)) {  // Optional, this check is also done by `if (missingValueReasons != 0)`.
    int missingValueReasons = 0;
    if (isNaN(v00)) missingValueReasons |= floatToRawIntBits(v00);
    if (isNaN(v01)) missingValueReasons |= floatToRawIntBits(v01);
    if (isNaN(v10)) missingValueReasons |= floatToRawIntBits(v10);
    if (isNaN(v11)) missingValueReasons |= floatToRawIntBits(v11);
    if (missingValueReasons != 0) {
        // At least one value is NaN.
    }
}

Listing F.8 — Example of code packing multiple reasons in a single NaN value

With this approach, declaring that a value is missing for multiple reasons is possible. For example, this is because of LAND and NO_PASS. A use case for this feature is: Consider a calculation involving many steps (resampling followed by convolution, etc.), each step using many inputs. If each step encodes for each output value all the reasons why the value is missing, and if a too large number of values appear to be missing at the end of the process chain, then users can check what is the main reason why values are missing by checking which NaN bit is most frequently set. For example, are final results missing mostly because of clouds? Or mostly because the remote sensor did not pass over the area? In the latter case, using a different remote sensing product may resolve the problem, while in the former case it may not unless the acquisition times are different enough.

F.3.4.  Code patterns leveraging NaNs

When codes can rely on NaN values behaving as specified by the IEEE 754 standard, some coding habits may be worth changing. Developers should keep in mind that any comparison ( <, >, <=, >= and ==) involving at least one NaN value will always be evaluated to false. This include x == x, which is false if x is a NaN value. Therefore, the following code:

void compare(double x, double y) {
    if (x != -9999 && y != -9999 && x < y) {
        // Do something if x<y AND both x and y are real values (not NaN).
    }
    if (x == -9999 || y == -9999 || x < y) {
        // Do something if x<y OR at least one of x and y is NaN.
    }
}

Listing F.9 — Comparisons using a Sentinel’s "no data" value

Can be replaced by the following code. Because the use of !(x >= y) may seem an unnecessary complication compared to (x < y) for developers who are not familiar with IEEE 754, a good precaution is to add a comment saying that the use of the negative form is intentional.

void compare(double x, double y) {
    if (x < y) {
        // Do something if x<y AND both x and y are real values (not NaN).
    }
    if (!(x >= y)) {    // Use `!` for entering in the block when a value is NaN.
        // Do something if x<y OR at least one of x and y is NaN.
    }
}

Listing F.10 — Comparisons relying on IEEE 754 behavior

F.4.1.  Getting the bits of a float

The following function is equivalent to java.lang.Float.floatToRawIntBits(float). Note that NaN values have a sign bit (it is possible to have “negative” NaN values), so codes really need the signed type shown below. The choice of signed or unsigned type changes the way that the max function will behave when checking which NaN has precedence in this test.

inline int32_t floatToRawIntBits(float value) {
    return std::bit_cast<int32_t>(value);
}

Listing F.13 — Float to raw int bits using the C++ 20 standard

Alternatively, bit patterns can be read directly from the data array without transiting by a floating-point type. This approach ensures that no Floating Point Unit (FPU) operation is involved, not even a copy (the truly paranoiacs could go further by applying reinterpret_cast on the array of bytes read from the file instead of the array of float values):

void foo(float* rasterData, int x, int y) {
    int32_t* dataAsInts = reinterpret_cast<int32_t*>(rasterData);
    int      offset     = y * rasterWidth + x;
    int32_t  bitPatern  = dataAsInts[offset];
    float    value      = rasterData[offset];
    // Do some stuff.
}

Listing F.14 — Float to raw int bits for a whole array

F.4.2.  Detecting NaNs without floating-point support

If the std::isnan(float) function is not available (for example, because of compiler options discussed in Annex F.4.4), quiet NaNs can still be detected from an array of IEEE 754 values interpreted as integers (e.g. using reinterpret_cast<int32_t*> as in Annex F.4.1) as below:

#define FIRST_QUIET_NAN 0x7FC00000

void foo(int32_t* dataAsInts) {
    int offset = ...;
    if ((dataAsInts[offset] & FIRST_QUIET_NAN) == FIRST_QUIET_NAN) {
        // The value is a quiet NaN.
    }
}

Listing F.15 — NaN check without floating point operation

For detecting both quiet and signaling NaNs, the mask would be slightly different, and another condition would need to be added for differentiating NaN from infinite values.

F.4.3.  Stability of standard functions

Some functions provided by standard libraries may have undefined behavior with NaN inputs. For example, the result of std::sort() on a vector<float> is only partially sorted when the vector contains NaN values. This issue is language dependent. For example, Arrays.sort(float[]) is fully defined in Java (NaN values are sorted last and −0 is sorted before +0), while Python behaves like C/C++. For languages with undefined behavior, users may need to filter the NaN values before invoking the standard function. Alternatively, users can sometime provide a custom comparison function. However, the filtering step or custom comparator is often needed anyway even with Sentinel’s “no data” values. For example, both NaNs and Sentinel’s “no data” values may need to be ignored when searching for the minimum or maximum values in an array because:

It may be unnecessary to filter the Sentinel’s “no data” values before sorting if the caller knows that these Sentinel’s values are smaller or greater than all valid values, but this is a special case. Filtering is also unnecessary in languages where the sort function is fully defined by default, including the NaN case.

F.4.4.  Effect of compiler options

By default, gcc and probably most C/C++ compilers produce code compliant with IEEE / ISO rules. The code generated with these default options passes the C/C++ test described in this section. However, the gcc compiler provides a -ffast-math option with the following warning:

The -ffast-math compiler option is a shortcut for the options listed in the following table. For verifying the effect of those options, the demo program has been compiled multiple times with the options added in the order listed below. The column on the right side reports whether adding an option made a difference compared with the previous compilation.

Table F.1 — Effect of -ffast-math compiler options

The tests pass with identical results for all options except -ffinite-math-only. With the latter option, the only behavioral change that impacted the demo program is in the std::isnan function. The sin, cos, tan, asin, acos, atan, atan2, floor, ceil, trunc, pow, sqrt, hypot and nanf functions continue to work as expected with NaN values. For working around the non-availability of the std::isnan function when the -ffinite-math-only option is used, the following check:

void foo(float value) {
    if (std::isnan(value)) {
        // Do some stuff.
    }
}

Listing F.16 — NaN check using std::isnan(…)

can be trivially replaced by the following, providing that missingValueReason is a signed integer and that all used NaN values are quiet NaNs and are “positive” (sign bit unset):

void foo(float value) {
    int missingValueReason = floatToRawIntBits(value);
    if (missingValueReason >= FIRST_QUIET_NAN) {
        // Do some stuff.
    }
}

Listing F.17 — NaN check without std::isnan(…) (limited applicability)

The above trick is possible in the test case (see the C/C++ source code) because of the way that missingValueReason is computed. This test uses max operations, together with the fact that “positive” quiet NaNs are greater than all other IEEE 754 values when compared as signed integers. Annex F.4.2 provides another way to identify the NaN values.

int main() {
    double x = std::sin(rand());    // The tested compiler options have no effect here.
    double y = std::cos(rand());
    double z = x + y + 0.0;         // Compilation result depends on -fno-signed-zeros.
    double w = z * y * 1.0;         // Multiplication by 1 is removed in all tested cases.
    printf("%f %f\n", z, w);
    if (z < w) {                    // Tests show no effect, unless < is replaced by ==.
        printf("z < w\n");
    }
    return 0;
}

ucomisd	(%rdx), %xmm0   # Accept quiet NaNs ('u' stands for "unordered").
jp	.L340           # Jump if unordered (i.e. if any operand is NaN).
jne	.L340           # Jump short if not equal.

comisd	(%rdx), %xmm0   # Signal an invalid operation for any kind of NaN.
jne	.L340           # Jump short if not equal.

F.5.1.  Expected results

The tests implemented in Java and C/C++ produce the same result, shown below. Those implementations use Math.fma(…) and std::fma respectively for “fused multiply-add” in bilinear interpolations. The drifts observed during the last iterations are caused by the intentionally chaotic nature of the test, which makes the last iterations very sensitive to small errors in previous iterations. This sensitivity is used for observing how large those drifts become when implementation details are changed or when C/C++ compiler options are added.

Errors in the use of raster data with NaN values in big endian byte order:
   Count     Minimum     Average     Maximum   Number of "missing value" mismatches
   11732      0,0000      0,0000      0,0000      0
   10771      0,0000      0,0000      0,0000      0
   10809      0,0000      0,0000      0,0000      0
   10762      0,0000      0,0000      0,0000      0
   10851      0,0000      0,0000      0,0000      0
   10814      0,0000      0,0000      0,0016      0
   10826      0,0000      0,0000      0,0871      0
   10918      0,0000      0,0042     18,1068      0
   10734      0,0000      0,0208     65,3695      9
   10792      0,0000      0,1692    127,6496     35

Listing F.24 — Drifts of test using double precision and FMA

The test implemented in Python differs from the Java and C/C++ implementations in that it does not use the FMA operation. Instead, the bilinear interpolation formulas are implemented in a more traditional way (except for + 0.0) as below:

v0 = (v01 - (v00 + 0.0)) * xf + v00
v1 = (v11 - (v10 + 0.0)) * xf + v10

Listing F.25 — Fragment of interpolation code in double precision (Python)

The reason for + 0.0 is because v00, v01, v10 and v11 are single precision floating-point numbers. The single precision nature of those numbers is explicit in Java and C/C++ code, but implicit and hard to see in the Python code. If the values are not cast to double, then (v01 - v00) and (v11 - v10) are computed with single precision. The precision lost causes a faster drift from the expected results. To force a cast to double precision, the Java and C/C++ code use the following:

double v0 = (v01 - (double) v00) * xf + v00;
double v1 = (v11 - (double) v10) * xf + v10;

Listing F.26 — Fragment of interpolation code in double precision (Java and C/C++)

In Python, adding 0.0 is a trick that appears to work. The result of the test in Python is shown below. The tests in Java and C/C++ produces an identical result if their code using FMA is replaced by the above code. This result shows a small loss of accuracy compared to the code using FMA. Note that on the machine used for the test (Intel® Core™ i7-8750H) and the gcc compiler used (14.2.1 20240912 (Red Hat 14.2.1-3)), the output is the same even with the -ffast-math option.

Errors in the use of raster data with NaN values in big endian byte order:
   Count     Minimum     Average     Maximum   Number of "missing value" mismatches
   11732      0,0000      0,0000      0,0000      0
   10771      0,0000      0,0000      0,0000      0
   10809      0,0000      0,0000      0,0000      0
   10762      0,0000      0,0000      0,0000      0
   10851      0,0000      0,0000      0,0000      0
   10814      0,0000      0,0000      0,0009      0
   10826      0,0000      0,0001      0,0871      0
   10918      0,0000      0,0047     18,1068      0
   10735      0,0000      0,0383     94,2899     10
   10789      0,0000      0,2025    127,6496     42

Listing F.27 — Drifts of test using double precision

If the single precision numbers are not forced to double precision before calculation, i.e. if the Python code is as below (as above but with + 0.0 removed):

v0 = (v01 - v00) * xf + v00
v1 = (v11 - v10) * xf + v10

Listing F.28 — Fragment of interpolation code in single precision

Then a larger drift is observed:

Errors in the use of raster data with NaN values in big endian byte order:
   Count     Minimum     Average     Maximum   Number of "missing value" mismatches
   11732      0,0000      0,0000      0,0000      0
   10771      0,0000      0,0000      0,0011      0
   10809      0,0000      0,0006      0,1469      0
   10760      0,0000      0,0264     12,2387      6
   10814      0,0000      0,3692    130,8832     90
   10670      0,0000      1,1468    138,7692    379
   10522      0,0000      2,2290    144,7353    799
   10379      0,0000      3,5549    158,2099   1383
    9938      0,0000      5,3378    161,3066   2001
    9783      0,0000      7,2493    171,4665   2595

Listing F.29 — Drifts of test using single precision

F.5.2.  Notes on an optimization strategy (optional)

TestNaN and TestNodata both use the same optimization strategy for selecting a missing reason in NO_PASS, LAND, CLOUD and UNKNOWN precedence order. This demo exploits the facts that:

The combination of these two facts allows the code to simply check for the maximum value, no matter if the raster has a mix of “no data” and real values. However, this trick would not work anymore if the code did not know in advance that all “no data” values are greater than all of the valid values. If they were less than the valid values, some >= operators would need to be replaced by <= in the TestNodata code (in addition of the max function being not applicable anymore). TestNodata would be even more complex if the raster had a mix of “no data” values that were both less than and greater than the real values (the use of abs could reduce this complexity but would require positive and negative “no data” values in the same range as the absolute values). By contrast, optimizations can be done more easily with NaN because the condition equivalent to #1 is more reliable.

G.2.1.  Experiment Datasets

Two types of data were selected for experiments of COG implementations. They are:

These datasets were selected for their relevance and utility in common geospatial applications, particularly in atmospheric and land surface studies. The objective of converting these datasets into Cloud Optimized GeoTIFF (COG) format is to enhance their accessibility and performance in cloud-based environments, enabling more efficient data analysis and application. In OGC Testbed 20, various schemes for implementing COG were experimented with to benchmark these typical data types in geospatial applications.

G.2.2.  COG Experiments

To optimize the conversion of Earth observation datasets into Cloud Optimized GeoTIFF (COG) format, various experimental parameters were tested. These parameters were chosen to evaluate their impact on performance, storage efficiency, and usability in cloud-based environments.

Block Sizes Tested:

Compression Methods: The primary focus was on LZW (Lempel-Ziv-Welch) compression due to its advantages in geospatial data applications:

G.2.3.  Implementation Details

The conversion of Earth observation datasets into Cloud Optimized GeoTIFF (COG) format was carried out with careful consideration of several key implementation details to ensure consistency, accuracy, and efficiency.

None selected

Skip to content
Using Gmail with screen readers
Enable desktop notifications for Gmail.
   OK  No thanks

Conversations
5.98 GB of 15 GB used
Terms · Privacy · Program Policies
Last account activity: 7 hours ago
Details
import os
from osgeo import gdal

def convert_to_cog(input_path, output_path):
    dataset = gdal.Open(input_path)
    if not dataset:
        print(f"Failed to open dataset: {input_path}")
        return

    translate_options = gdal.TranslateOptions(
        format='COG',
        creationOptions=[
            "COMPRESS=LZW",
            "BLOCKSIZE=512"
        ]
    )
    gdal.Translate(output_path, dataset, options=translate_options)
    dataset = None  # Close the dataset

def convert_nc4_to_cog(input_file, subdatasets, base_output_dir, date_str):
    output_dir = os.path.join(base_output_dir, "MERRA2", "MERRA2-512M-LZW")
    os.makedirs(output_dir, exist_ok=True)  # Ensure output directory exists

    for subdataset in subdatasets:
        subdataset_path = f'HDF5:"{input_file}":{subdataset}'
        output_file = generate_nc4_output_filename(output_dir, subdataset, date_str)
        try:
            convert_to_cog(subdataset_path, output_file)
            print(f"Converted {subdataset} to COG: {output_file}")
        except Exception as e:
            print(f"Failed to convert {subdataset}: {e}")

def generate_nc4_output_filename(base_dir, subdataset, date_str):
    variable_name = subdataset.strip("//")
    filename = f"MERRA2_{variable_name}_{date_str}_512m_COG.tif"
    return os.path.join(base_dir, filename)

def generate_tif_output_filename(base_dir, filename):
    base_name = os.path.splitext(filename)[0]
    output_filename = f"{base_name}_512m_LZW_COG.tif"
    return os.path.join(base_dir, output_filename)

def process_tiff_files(input_folder, base_output_dir):
    output_dir = os.path.join(base_output_dir, "HLS.L30", "HLSL30-512M-LZW")
    os.makedirs(output_dir, exist_ok=True)  # Ensure output directory exists

    print(f"Processing TIFF files in folder: {input_folder}")
    for filename in os.listdir(input_folder):
        if filename.endswith(".tif") and not filename.startswith("._"):
            input_path = os.path.join(input_folder, filename)
            output_path = generate_tif_output_filename(output_dir, filename)
            print(f"Processing TIFF: {filename}")
            try:
                convert_to_cog(input_path, output_path)
                print(f"Converted {filename} to COG: {output_path}")
            except Exception as e:
                print(f"Failed to convert {filename}: {e}")

if __name__ == "__main__":
    # Input paths
    nc4_file = "/Volumes/T7 Shield/CSISS/OGC_Testbed_20/Data/MERRA2_400.tavg1_2d_slv_Nx.20140101.nc4"
    input_folder = "/Volumes/T7 Shield/CSISS/OGC_Testbed_20/Data/Landsat/HLSL30_2.0-20240912_144107"
    base_output_dir = "/Volumes/T7 Shield/CSISS/OGC_Testbed_20/Data/COG_converted"

    # Ensure base output directory exists
    if not os.path.exists(base_output_dir):
        os.makedirs(base_output_dir)

    # Subdatasets for NC4 files
    subdatasets = [
        "//CLDPRS", "//CLDTMP", "//DISPH", "//H1000", "//H250", "//H500", "//H850",
        "//OMEGA500", "//PBLTOP", "//PS", "//Q250", "//Q500", "//Q850", "//QV10M",
        "//QV2M", "//SLP", "//T10M", "//T250", "//T2M", "//T2MDEW", "//T2MWET",
        "//T500", "//T850", "//TO3", "//TOX", "//TQI", "//TQL", "//TQV", "//TROPPB",
        "//TROPPT", "//TROPPV", "//TROPQ", "//TROPT", "//TS", "//U10M", "//U250",
        "//U2M", "//U500", "//U50M", "//U850", "//V10M", "//V250", "//V2M", "//V500",
        "//V50M", "//V850", "//ZLCL"
    ]

    # Process the NC4 file
    print(f"Processing NC4 file: {nc4_file}")
    date_str = nc4_file.split('.')[-2]  # Extract date from filename
    convert_nc4_to_cog(nc4_file, subdatasets, base_output_dir, date_str)

    # Process all TIFF files in the input folder
    process_tiff_files(input_folder, base_output_dir)

G.3.1.  Resulted COG Datasets

The resulting COG formatted datasets are shown in the following table — Table G.1. Two types of datasets were experimented with, including one from the MERRA-2 dataset and another from the Landsat 9 HLS sensor. Three types of block sizes were applied to both datasets: 128 by 128, 256 by 256, and 512 by 512.

Table G.1 — Output images comparison

	Merra2	HLSL-30
Original
128 x 128
256 x 256
512 x 512

G.3.2.  File Sizes

Table G.2 — Sizes of File

G.3.3.  Read Experiments

The following table (Table G.3) shows the results of reading experiments.

Table G.3 — Reading of COG Files

Observations from the Reading Experiment:

G.3.4.  Tiling and HTTP Range Request Support

The tiling and HTTP range request capabilities of the COG files were inspected using gdalinfo. Below is an example output for the 256×256 block size configuration:

gdalinfo MERRA2_example.tif | grep "Block"
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray

Listing G.2

Metadata Summary:

H.3.1.  Analysis

In August 2024, a short, geotagged video was filmed aboard a yacht under sail using an Android device running Away Team’s TrkdCam application. The data were automatically timestamped using the WebVMT media start time in Coordinated Universal Time (UTC). Previous synchronization work in OGC Testbed 19 (OGC 23-042) quantified typical timing errors to be within 2 seconds due to device latencies.

H.3.1.1.  Data Aggregation

AIS data were downloaded from a vessel tracking website using the Maritime Mobile Service Identity (MMSI) of the yacht. These data included the vessel’s location, heading and speed at intervals varying between 3 and 10 minutes over a period of around 6 hours. Satellite-AIS data were not available, so the yacht’s trajectory was divided into two segments that correspond to the outward and return legs of the trip.

Open source data from OpenSeaMap were integrated with the WebVMT JavaScript playback engine to display nautical chart data including seamarks.

A new utility was developed from the GPX Importer community tool to import AIS data in a Comma Separated Value (CSV) file into WebVMT format. WebVMT code developed in OGC Testbeds 17 and 19 was used to synchronize and aggregate imported AIS data with on-board video metadata using the clip and merge functions.

 

Figure H.1 — AIS Data Synchronized and Aggregated With On-Board Video

The AIS trajectory in Figure H.1 is shown as a faint blue line for the outward journey and a stronger blue line for the return trip. AIS data locations are marked along the line to highlight their sparse temporal and spatial nature. On-board data is shown in orange in Figure H.2.

 

Figure H.2 — Comparison of AIS Trajectory (Blue) With On-Board Trajectory (Orange)

The location of the Lepe Spit mark — shown in red in Figure H.2 — was added from the 2024 GPX data published on the Solent Cruising and Racing Association (SCRA) website and circular zones were added with 20m and 50m radii for scale. Circles were also added to the AIS data location with 10m and 50m radii to help gauge distance.

There is a discrepancy between the mark’s 2024 location and its OpenSeaMap position, though some margin of error should be allowed for tidal effects and the locations are still within 20 meters. AIS and on-board trajectories are within 10 meters once they pass the AIS data point in Figure H.2. The tracks diverge prior to this data point so further analysis was required.

H.3.1.2.  GEOINT Analysis

At the start of the footage, the Lepe Spit mark is less than 100m away and is clearly visible in the video imagery. This mark closely resembles the symbol indicating its location on the nautical chart. Information published on the UK Trinity House website confirmed that this is a South Cardinal mark designed to warn of a maritime hazard to its north, and advises that “mariners will be safe if they pass…​south of a south mark.” An area of shallows is marked to the north of its location on the OpenSeaMap chart and therefore it is reasonable to assume that the vessel passed to the south of this mark.

The interpolated AIS track passes north of this mark in Figure H.2, so a single additional data point was added to the trajectory to correct this error. This location was estimated to be 20m from the Lepe Spit mark location to allow a margin of error for safe navigation. The trajectory time was estimated from linear interpolation of the two adjacent AIS data points assuming movement over ground as a first approximation.

 

Figure H.3 — Comparison of Corrected AIS Trajectory (Blue) With On-Board Trajectory (Orange)

The corrected AIS trajectory shown in blue in Figure H.3 closely matches the on-board data shown in orange — both in terms of location and time. This estimate could be refined to model contributions from tide, engine, wind and sail, but the simple approximation used here is consistently accurate to within 10m and serves to demonstrate how basic aggregation can significantly enrich geospatial data.

H.3.2.  Results

Geospatial data from multiple sources were aggregated and accurately synchronized with video using WebVMT. These data were successfully integrated with other web technologies in a web browser to help guide GEOINT analysis and were significantly enriched using simple techniques. Metrics show that the results obtained are accurate to within ten meters and synchronized to within two seconds.

Figure H.4 — Corrected AIS Data (Blue) Aggregated With On-Board Data (Orange) Synchronized With Video

WebVMT tools developed in OGC Testbeds 17 and 19 were reused and updated to support negative cue times for non-destructive synchronization. A new utility was created to import GPX data — based on the online community tool for importing GPX data — and extended to read CSV files to facilitate data migration to WebVMT.

OpenLayers API was integrated with the WebVMT JavaScript playback engine to support OpenSeaMap nautical chart data. A new engine feature was prototyped to generate an event after the map element has been initialized in the browser which enables the OpenSeaMap seamark layer to be properly integrated.

The prototype code utilizes the HTML Event interface (CustomEvent) to trigger execution of user code, instead of allowing that code to be called directly by the current thread. This design enables proper sandbox isolation that mitigates any threat from malicious code in the webpage and is consistent with HTML security guidelines.

H.4.1.  Analysis

A video file with in-band metadata encoded in MISB ST 0601 was identified as a good candidate with which to investigate camera orientation issues. The footage was filmed from an aircraft that was circling over Wyoming, USA with a video camera trained on a truck that was travelling along a highway below. The relative motions of the aircraft and truck result in a complex sequence of camera angles that are constantly changing.

A C++ utility developed in OGC Testbed 16 was modified to export the SMPTE Key Length Value (KLV) sensor location (MISB ST 0601 tags 13, 14 & 15) and frame center (MISB ST 0601 tags 23, 24 & 25) video metadata to generate two WebVMT paths. These trajectories represent the tracks of the camera and its subject location respectively.

H.4.1.1.  Camera Heading

Map rotation support is included in several publicly accessible Web Map APIs including OpenLayers and MapLibre. The WebVMT JavaScript playback engine was updated to include the MapLibre API in order to prototype this feature.

 

Figure H.5 — OpenLayers Map Rotation Aligned With Video Camera Heading

A WebVMT heading attribute was prototyped in OGC Testbed 19 for the sight lines used to track vehicles from video footage (OGC 23-042). These data were also included in the Hillyfields Bubble dataset that was published under an open source license in OGC Testbed 19, and subsequently proposed as a new feature in the WebVMT Community Group (CG). The CG proposal was extended to include a heading attribute for map rotation and implemented as a prototype feature in the WebVMT engine.

A new utility was developed to calculate the heading of one WebVMT path from another. This was used to create WebVMT content centered on the truck trajectory with the map rotation automatically aligned to the camera heading.

 

Figure H.6 — MapLibre View Aligned With Video Camera Orientation

Figure H.6 shows the frame center track on the map in blue which is marked by a white crosshair baked into the video frame above. This is a proxy for the truck’s trajectory in this case.

 

Figure H.7 — Comparison Between Paths For Camera (Orange) & Frame Center (Blue)

There are two ways in which to verify that the map rotation accurately reflects the camera orientation in this example:

1. The video includes a compass indicator which is baked into the image on the left-hand edge towards the top of the video frame in Figure H.7. This indicator matches the compass direction shown in the top right corner of the map element. These indicators are also visible in Figure H.6 and Figure H.5.

2. Figure H.7 shows the paths of the camera in orange and frame center in blue. The tip of the blue path is always vertically above the tip of the orange path which also confirms that the map rotation is correct.

H.4.1.2.  Cognitive Analysis

Loading these same video and VMT files into a web page with photogrammetric data can assist cognitive analysis of features that are visible in the video, though may not be marked on a map.

 

Figure H.8 — Truck Video With Rotated Photogrammetry Data

The top video panel in Figure H.8 shows the truck passing under a sign that spans the road on a gantry. The bottom panel shows the photogrammetry data synchronized to the video with WebVMT which includes several notable features:

1. The sign gantry is visible in the center of the photogrammetry data and can be recognized from the shadow that it casts. This prominent feature provides a simple check that the data are probably in the right location.

2. A paler area of ground that is parallel to the highway and next to the gantry can be seen beyond the fence line in the foreground of the video frame. This paler area extends to the near side of the fence line to form a shape that is similar to the letter ‘J’. A feature with this same unique shape can be seen in the lower map panel below the end of the blue line which provides good evidence that these two images show exactly the same location.

3. The tip of the frame center trajectory shown in blue in the map panel accurately matches the location of the crosshair in the video frame above. The relative location of the trajectory tip to the foot of the gantry support indicates that this is accurate to within a meter or so. This provides a metric of the data accuracy and demonstrates good synchronization.

4. Having established that the video and photogrammetry images show the same location, less obvious features can also be identified. Crash barriers can be seen on either side of the near carriageway which lead up the sign gantry and then end a few meters beyond this point. These features are visible in both images and their locations can be accurately identified.

5. The inside lane markings change from short dashed lines to a continuous white line on both carriageways shown on the right side of the video image. The same changes in lane markings can be seen towards the right of the photogrammetry data. The locations of these points can be accurately identified.

H.4.2.  Results

Aligning map rotation with camera orientation improves geospatial cognition of geotagged images and simplifies GEOINT analysis. Vehicle satnav systems and smartphone navigation applications have adopted the same proven technique to display geospatial information in a concise form that is aligned with the user’s heading.

A new utility was developed to calculate the heading of one WebVMT path from another to support camera orientation. The data produced by this tool successfully demonstrated that WebVMT can synchronize and interpolate dynamic data, including camera orientation, to match video imagery with frame-level accuracy. Accurately synchronizing photogrammetry data with a video frame enables location verification, provides accuracy metrics and facilitates cognition of details that might otherwise have been overlooked in GEOINT analysis.

H.5.1.  Data Capture

In April 2024, Away Team collaborated with Ordnance Survey (OS) at their headquarters site near Southampton, UK to capture geotagged video for a roadside litter monitoring use case.

Target litter items were placed at known roadside locations on site and filmed using vehicles equipped with dash cams. The litter included tins, drink cans and plastic bottles in a variety of sizes and colors, and traffic cones. The cones were identified in OGC Testbed 19 (section 7.4, OGC 23-042) by the UK National Highways Agency as an item of interest for collection.

Two vehicles were equipped with front- and rear-facing cameras to capture geotagged video:

Figure H.9 — Ordnance Survey StreetDrone Equipped With Extra Cameras

Five test runs were completed to capture data from the same road circuit. Alternate vehicles were used for each successive run.

H.5.2.  Analysis

Each vehicle produced two separate geotagged video files corresponding to its front and rear cameras. Video output from the Nextbase system was automatically segmented into files with a maximum duration of one minute. This automatic feature allows short video files to be submitted as police evidence without any editing process to ensure that their forensic integrity is preserved.

Figure H.10 — Nextbase Front And Rear Cameras In Situ

The first challenge was to handle the file segmentation process seamlessly.

H.5.2.2.  Litter Categorization

The next challenge was to identify litter observations from the video footage.

Classification of litter items is logistically important to ensure that the assigned collection team is properly equipped to remove the identified items. Glass bottles are heavier and bulkier to carry than their plastic counterparts, and a broken glass bottle is more hazardous to collect than a plastic one.

 

Figure H.11 — Target Litter Items Used In Data Capture

Reference photos were used to create a catalogue of litter locations and classify each item with a unique identifier such as ‘clear-plastic-bottle-01’. This allowed the same target item to be identified from multiple observations.

Being able to categorize and count litter items is important to successfully monitor accretion and deploy collection teams efficiently. Tracking individual small items such as bottles or cans is impractical and unnecessary. However, tracking larger objects such as a mattress or a wooden pallet is important since these items can pose a significant threat to the safety of road users. Larger items may require specific equipment to safely collect and remove them, so this information is vital for collection team deployment.

H.5.2.5.  Data Synchronization

WebVMT negative cue times were prototyped in OGC Testbed 19 (section 7.4, OGC 23-042) for non-destructive synchronization of data. These were used in OGC Testbed 20 to synchronize the cumulative dash cam trajectory data with each video segment file.

 

Figure H.12 — Front Video Segment With Litter Observations (Red) & Dash Cam Trajectory (Blue)

Figure H.12 shows a video segment synchronized with the camera trajectory — in blue — which extends prior to the start of that file without loss of data integrity.

Litter observations were also synchronized and aggregated with the camera trajectory. An active observation at the start of a new video segment file can be recorded with a negative cue start time. This prevents truncation due to segmentation and ensures that data integrity is retained.

 

Figure H.13 — Rear Video Segment With Litter Observations (Red) & Dash Cam Trajectory (Blue)

Figure H.13 shows a target litter item — in red — observed by the rear camera. This same item was also observed by the front camera in Figure H.12.

These data were ready for further study but insufficient time remained in OGC Testbed 20. This analysis could be deferred to OGC Testbed 21 as described in Future Work.

H.5.3.  Results

A number of video players were evaluated for use as an inspection tool and compared against the requirements identified for metadata analysis. These features are:

Table H.1 — Key Analysis Features Supported by Video Players.

The mpv application offered the best support for the identified analysis requirements. Zoom and pan support were provided by adding the following customized keyboard shortcuts:

q   add video-zoom -.25
e   add video-zoom .25

w   add video-pan-y .05
d   add video-pan-x -.05
x   add video-pan-y -.05
a   add video-pan-x .05

s   set video-pan-x 0; set video-pan-y 0; set video-zoom 0

Listing H.1 — Zoom/Pan Keyboard Shortcuts For mpv (input.conf)

The IINA application is based on mpv and offered comparable support for the analysis requirements but was hindered by a bug in the timing display.

Factors affecting observation visibility in the video were identified including field of view, obstruction by foreground and background objects, and video resolution and compression.

Front and rear litter observations have been geotagged in segmented video files without loss of data integrity. This dataset is ready for further study but insufficient time remained in OGC Testbed 20 to complete this analysis. This analysis could be deferred to OGC Testbed 21 as described in Future Work.

H.6.1.  GIMI Design Considerations For Video Metadata

In-band and out-of-band metadata encodings are different approaches to the same problem that are not mutually exclusive. Both have strengths and weaknesses, and the best choice should be determined from analysis of the individual use case.

NOTE KLVs encoded in WebVMT

00:00:12.345 -->
{ "sync": { "type": "org.ogc.testbed-20.gimi.example.klvs", "data":
  { "key-1": "value-1",
    "key-2": "value-2",
    "key-3": {"key-3a": "value-3a", "key-3b": "value-3b" } }
} }

WEBVMT

MEDIA
url:my-video.mp4
mime-type:video/mp4
start-time:2018-02-19T12:34:56.789Z

NOTE Zero on the media timeline for the my-video file corresponds to UTC timestamp 2018-02-19T12:34:56.789Z

00:00:00.000 -->
...

H.6.2.  Video Search Benchmark Analysis

WebVMT offers a common format for sharing timed video metadata on the web that is agnostic of video encoding and accessible to web crawlers. This enables online search engines to handle metadata-related queries for video data such as GIMI files.

Efficient access to metadata-related content in a video file is an important benchmark for the proposed GIMI design and serves to highlight key differences between the in-band and out-of-band metadata approaches.

The search process has two stages:

These two steps may be combined for in-band metadata use cases to help mitigate the access overhead. Selective query benchmarks may have similar file access requirements to the first benchmark. However, they also serve to quantify the relative sizes of the metadata access overhead, query selection and retrieval steps in the overall process.

To compare in-band and out-of-band queries, identical metadata should be encoded into both formats for direct comparison.

H.6.3.  Results

Use of out-of-band video metadata can improve data privacy and security since the sidecar file has security permissions that are separate from the video file. This enables the two files to have different access permissions. More open permissions may be granted to access the sidecar file which enables search queries to determine whether the video file is of interest. This approach does not require direct access to the video content which may contain personally identifiable information such as images of faces or vehicle registrations.

Out-of-band metadata search queries only require access to the associated video file when a positive result is returned. Monitoring this file access pattern could identify malign access and help to curb web scraping activities employed by unscrupulous organizations abusing Big Data and Artificial Intelligence (AI) technologies.

In addition, use of out-of-band metadata queries should reduce file access bandwidth since a negative result only requires access to the smaller sidecar file and should never access the larger video file.

No benchmark metrics were calculated in OGC Testbed 20 due to the lack of suitable GIMI video files for analysis.

H.7.1.  Maritime

Away Team has successfully demonstrated how geospatial data from multiple sources can be aggregated and accurately synchronized with video using WebVMT. Web browser integration with web technologies such as OpenSeaMap has made data more accessible, helped to guide GEOINT analysis and significantly enriched data using geospatial techniques.

H.7.2.  Camera Orientation

Combining camera orientation information with map and photogrammetry data can significantly improve cognition of imagery including features that are not normally marked on maps such as trees and road furniture. This demonstration highlights the importance of camera orientation data in the GIMI standard for still and motion imagery, and how GEOINT analysis can benefit from this information.

H.7.3.  Litter Monitoring

Video player requirements have been identified for geotagged video analysis including millisecond display, frame-level control, playlist support for segmented video, and zoom and pan to identify visual details within video frames.

Factors affecting observation visibility in video have been identified including field of view, obstruction by foreground and background objects, video resolution and compression.

The litter monitoring dataset is ready for investigation but there was insufficient time in OGC Testbed 20 for further study. This analysis could be deferred to OGC Testbed 21 as described in Future Work.

H.7.4.  GIMI Video Metadata Considerations

MISB Class 3 Imagery is an increasingly important resource due to proliferation of location-aware devices with cameras such as drones, dash cams, body-worn video and smartphones.

In-band and out-of-band metadata have strengths and weaknesses which depend on the individual use case. These two approaches are not mutually exclusive, and GIMI video metadata design should address both options.

JSON encoding of SMPTE KLVs can be facilitated by WebVMT which is well-integrated with HTML and CSS for access, sharing and searching online.

TAI and UTC timestamps can be integrated with GIMI for in-band and out-of-band video metadata.

H.7.5.  Video Metadata Search

Out-of-band metadata allows web crawlers to access video metadata on the web using sidecar files in a format that is agnostic of the video encoding. This enables online search engines to support video metadata queries.

Use of out-of-band metadata queries can improve privacy and security of personally identifiable information contained in video files, help to identify malign access and reduce file access bandwidth.

Benchmark tests have been designed to provide metrics for video metadata search access with identified use cases.

H.7.6.  WebVMT Enhancements

Analysis of these video use cases has highlighted new requirements for the WebVMT design which have been implemented as prototypes in the engine and tools.

OpenLayers and MapLibre APIs have been integrated with the WebVMT JavaScript playback engine to support a wider range of map data and interface capabilities.

A new WebVMT feature has been prototyped to facilitate integration with maps that have customization requirements such as OpenSeaMap.

A new WebVMT utility has been created to calculate the heading of one path from another to support camera orientation.

WebVMT tools developed in OGC Testbeds 17, 19 and 20 have been updated into a single Perl command-line utility designed to create and manipulate WebVMT files to facilitate video metadata analysis.

The online community tool for importing GPX data has been ported to a command-line utility and extended to read CSV files to help simplify data migration to WebVMT.

Negative cue times and unbounded cues are breaking changes to WebVTT and the text/vtt MIME type. A new IANA media type text/vmt should be considered to support these HTML features and improve timed video metadata integration with web browsers.

H.8.1.  Data Analysis of Litter Monitoring Use Case

Analyze data gathered in OGC Testbed 20 to demonstrate how roadside litter can be monitored using dash cam footage from a vehicle fleet.

H.8.2.  Benchmark Metadata Search Use Cases With GIMI Video

Identify suitable files for video metadata benchmark testing and calculate metrics to compare in-band and out-of-band search request to: