desruisseaux
commented
3 years ago Closing this issue since I posted it on the wrong project. Sorry for the disturbance, and thanks to Roger for moving it to the right place.
Closed desruisseaux closed 3 years ago
desruisseaux
commented
3 years ago Closing this issue since I posted it on the wrong project. Sorry for the disturbance, and thanks to Roger for moving it to the right place.
CRS definitions
Source and target CRS are specified either by ISO 19162 WKT strings or by “authority:version:code” tuples (version is optional). GML may be an option too, but it is not up to date with latest ISO 19111.
CRSs do not need to be geographic; any ISO 19111 CRS is okay, including three- or four-dimensional CRS. Two dimensional geographic and projected CRS may be the minimum that all implementations shall support, while CRS with more than 2 dimensions may be considered optional features.
Source and target CRS shall have the same number of dimensions (or should they? It may be possible to relax this restriction, but that could be another discussion).
Coordinate system axes
Source and target CRS should have axes in the same order and same direction. We could also add a “same units” restriction, but the simplification brought by this restriction is small.
Alternative: we could relax above constraint by introducing interpolation CRS. I'm not sure if it would be exactly the same than ISO 19111 interpolation CRS, but I think it is close in its intent. The set of rules could become:
The intent is to allow some formulas (scale, rotation, logarithms, …) to be applied before to interpolate in the grid.
Grid shape
The grid shall be a rectangle (2D), a cube (3D) or an hyper-cube (4D+). The restriction to rectangular shape is for performance reasons. Testing whether a point is contained inside a shape would be much more expensive if arbitrary polygons were allowed.
Grid nodes shall be equally spaced. This constraint is mentioned here for clarity, but it is actually a mathematical consequence of the 2 next sections.
Grid coordinates
Each node is identified by integer coordinates. Grid node coordinates start at zero and increase by 1 for each node. The example below shows coordinates for a 4×4 grid:
Note: (0,0) is not necessarily in the geographic north-west corner of the grid — see next section.
Conversion from CRS coordinates to grid coordinates
We need to define either a “grid to CRS” or a “CRS to grid” conversion. The following discussion takes the “grid to CRS” convention by analogy with raster data such as TIFF “world files”, but switching between conventions is just a matter of inverting matrices.
The conversion from grid node coordinates to interpolation CRS (or source CRS) is defined by an affine transform (e.g. the “tfw” file associated to a TIFF “world file”). For example if grid node coordinates (0,0) correspond to geographic coordinates (λmin, φmin) and if each node is separated from next node by ∆λ degrees of longitude and ∆φ degrees of latitude, then conversion from grid node coordinates (i,j) to geographic coordinates (λ,φ) can be expressed as below:
λ = λmin + ∆λ ⋅ i φ = φmin + ∆φ ⋅ j
In matrix form (the last row is often omitted in operation parameters or in TFW files):
If the CRS has (latitude, longitude) axis order instead than (longitude, latitude):
If the grid (0,0) coordinates map to upper-left corner instead than lower-left corner:
Specifying affine transform coefficients allow to handle various axis orders and directions in non-ambiguous way. It is easily generalizable to 3D or 4D cases (we simply add rows and columns) and allows handling of rotations if desired. Graphical libraries work with affine transforms for decades and modern hardware such as GPU are designed for manipulating them efficiently.
If handling arbitrary affine transforms is considered too complex, we could define which kind of matrix shall be supported by all implementations (e.g. matrices having only scales and translations terms, without rotation) and said that supporting more general matrices is optional.
Domain of validity
A point is considered inside the grid if it is inside the smallest bounding box containing all grid nodes. So the valid bounding box is the yellow area on the left side, not the right side. The right side is mentioned because it is a common assumption when working with raster data, where values are often assumed located at pixel centres. In raster case the domain of validity is expanded by ½ cell size on each side. This expansion shall not be applied on GGXF grids.
Above definition assumes that geographic coordinates (φ,λ) are converted to grid coordinates (i,j) before to determine whether the point is inside the grid. This is necessary for handling the general case of affine transform, where a rotation may exist. However implementations are free to test directly on geographic coordinates if desired. It is easy to do if the “grid to CRS” affine transform contains only scale and translation terms (reminder: as proposed in previous section, we may say that only such transforms are mandatory in core standard).
Values at nodes
Each node contains an arbitrary amount of values. In a GeoTIFF file each value may be stored in a band. In a netCDF file each value may be stored in a variable. The set of “bands” or “variables” is called “sample dimensions” (not to be confused with CRS dimensions) in ISO 19115. In a GGFX file they could be:
Above is only an example. The set of legal sample dimensions is to be determined.
Coordinates in target CRS can be computed by adding an offset to coordinates in source (or interpolation) CRS. It is advantageous compared to storing target coordinates directly because achieving a centimetric precision on latitude and longitude values require the use of double-precision floating point values, while offsets require only single-precision floating point values if they are small displacements. It divides by 2 the grid size. In addition the use of offsets is more convenient for computing the inverse transform (from target CRS to source CRS).
Null values
Grid should contain a value at every nodes. The situation is different than raster data, where it is common to have missing values (e.g. Earth surface behind clouds, etc.). In case of raster data, null values are rendered by assigning a different color to the corresponding pixel. But in the case of GGXF files, null values can complexify the task of identifying which pixels are affected. It also make inverse transforms with iterative algorithm like NTv2 more hazardous if the point to transform is valid but close to a null value. When a shape crosses an area with null values, the software may need to connect the dots in some way. Consequently, implementations may be forced to replace null values by some approximated values before to do any work. Leaving that work to implementations creates a risk of discrepancies between different software using different replacement strategies, and the result would not be as reliable as if this work was done by the producer.
I propose to require that some reasonable value is provided in all cells, and flag “null” values with a high uncertainty. Maybe the metadata could contain a sentence like “every node with an uncertainty greater than 1 meter should be considered as a missing value”.
Uncertainty can be shown in various ways by software: the status bar showing mouse coordinates can contain a “± 5 cm” suffix. Or geometries may be drawn with a line of variable thickness, where the line thickness is the uncertainty at that location.
Interpolations
Conversion from geographic coordinates (φ,λ) to grid node coordinates (i,j) may result in fractional i and j coordinate values, in which case an interpolation is required. Bilinear interpolation is made simple by the fact that the interval between each grid node is exactly 1.
It is difficult to encode formulas in a data file. The interpolation method may need to be specified only by a name (e.g. “bilinear” or “bicubic”), in which case we may need to enumerate them in the specification (or enumerate at least what should be the core set).
Extrapolations
Implementation should not extrapolate values outside the grid. For grids of offsets value (e.g. NTv2), I suggest to take the offsets of the closest point in the grid. I suggest to not set the offset to zero (unless the closest point in the grid has zero offsets) for avoiding discontinuities when a line crosses the border between inside and outside the grid. Some high uncertainty (to be determined) would implicitly apply to all points outside.
Continuity at grid border
If two grids share a common border, the above “domain of validity” definition implies that both grids have nodes located on the same border. This is shown by the green dots in the picture on the right side. This is different than raster data shown on the left side, where having a common border does not imply that both rasters provide values at the same location. In other words, the GGXF grid convention implies nodes overlap (the green dots) that do not exist in the raster case.
However the GGXF convention avoids the difficulty of interpolating in the area which is half yellow and half blue in the raster case. This convenience come at a cost, which is that a discontinuity may exist if the “yellow nodes” and the “blue nodes” below green dots do not have the same values.
More generally, when interpolating values along a line from left to right, the interpolated values should not “jump” by an amount greater then uncertainty at the point where we move from yellow grid to blue grid.
Multiple grids
A file can have any number of base grids. We do not need to require the base grids to be non-overlapping, or to put constraint about their relative position. However doing so may enable some processing. For example an implementation may want to merge all base grids in a single huge virtual grid. I’m not sure that a lot of implementations would want to do that, but putting above-cited restrictions would enable this scenario. This is a case where putting more restrictions on one side gives more flexibility on the other side.
A base grid can have any number of child grids. Each child grid should be fully contained in its parent (for making easier to build R-Tree). If this is not possible, the “child” grid can be defined as a standalone (possibly overlapping) base grid, but producers are encouraged to keep the amount of base grids low because a tree-like structure like “parent-children” relationship can be traversed more efficiently.
Except for above-cited containment rule, the restrictions documented in NTv2 format are not strictly necessary. But enforcing them enable some scenarios similar to the way we build pyramid of images. A possible compromise may be to take NTv2 rules as recommendations instead than requirements.