GGXF grid structure proposal

The post below is from Martin Desruisseaux on https://github.com/Esri/ggxf/issues/3

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: • Interpolation CRS shall have axes in the same order and direction than target CRS. • No such constraint on source CRS, unless the interpolation CRS is the source CRS (probably the most common case). • If interpolation CRS is not the source CRS, then an ISO 19111 conversion from source CRS to interpolation CRS shall be provided.

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:

Offset to add to the first coordinate in source (or interpolation) CRS.
Offset to add to second coordinate.
… same for each additional dimension, if any …
Uncertainty at that node (unit specified in metadata). 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.

This looks a great start.

Null values.

I feel (as per deformation model discussions) that there is merit allowing for a null/undefined value. If determining a value at a point requires using an null value, then that value is null. This requires software to return some equivalent of undefined. Depending on the language this might NaN (not a number), null, etc. Alternatively the software could throw an exception. Similarly this could be the case beyond the extent of the grid. If we permit null values then the GGXF metadata should define the behaviour at points beyond the grid (null or zero).

Continuity

There is an implication that the spatial function implemented by the grid should be continuous (with the exception of null values if these are permitted). Within a grid the interploation functions we are likely to use will ensure that. However this is not guaranteed at grid boundaries where the value beyond the edge of the grid is not null (eg adjoining, overlapping, or nested grids). Do we specify that it must be for a valid GGXF grid, or do we allow the possibility of discontinuity and provide metadata to inform software if the grid is continuous.

Also the discussion above on continuity at a grid border proposes that the adjoinging grids can be treated independently. However if we are using eg bicubic interpolation then it would be more consistent to use nodes from both grids for the rows near the edge (or overlap a row/column at the edge)

Multiple grids A file can have any number of base grids. ...

I note t that this has only a start mutiple grid models. There are a number of factors to consider. For example:

what is a base grid
if there are overlapping base grids how should the correct value be defined. Is one selected, do we define a priority order for selecting potential grids, Do we assume they will evaluate to the same value (tricky if we are not using bilinear). Could/should the values from overlapping grids be added together rather than one selected,

Nested grids may have a number of options, eg

do they replace or add to the parent grid values where they overlap
must the values be continuous across the boundary the the parent grid (or possibly with an adjoining child grid, in which case it need not be continuous with the parent grid).
if the nested grid replaces the parent, and if we require continuity, then that puts some constraints on the grids. The child grid must subdivide the parent grid, not be offset from it. Bicubic interplation cannot be supported.
if the nested grid values add to the parent grid value these constraints are not necessary - all that is required is that the child grid evaluates to zero at its edge (except where it is adjacent to another child grid)
can siblings overlap, how should this be handled, which is similar to the question of overlapping base grids.

There are a lot of options!

GGXF needs to specify which of these are supported, what conditions a valid grid will meet, and what additional conditions it may offer to meet by metadata specifications. For example we could specify that a valid grid may be discontinuous but all grids have a flag specifying continuous or otherwise. Software could then implement a subset - eg only recognize continous grids for some applications.

Thanks Roger and Chris. There is plenty of food for thought and issues to resolve. It is also beneficial that discussions from the DMFM work is being cross-considered. I can only bit off a little at at time here, so starting at the beginning:

CRS definitions Two dimensional geographic and projected CRS may be the minimum that all implementations shall support...

The following may seem obvious, but I'll pose the question anyway. While affine transformation for projected CRS can be defined, they are not shown or mentioned in the section "Conversion from CRS coordinates to grid coordinates". Is it the intention that GGXF allow for any applicable coordinate operation to be performed directly in projected space, or is there an unstated assumption that conversion to geographic space should occur? Also, some statements should be made about the differences, if any, and methodologies of applying/using GGXF in projected space. This is merely a sanity check to ensure GGXF will be equally applicable in projected as well as in geographic space.

Do the rules under Coordinate system axes prevent or allow mixture of projected and geographic space for source, target, and interpolation CRS's? Not sure from the wording.

Hello Chris. Below is my 2 cents:

Null values

Returning null, NaN or throwing an exception is easy when transforming points without worrying about their connections with other points. But a difficulty happens when the point is part of a geometry. A null point may force to split a polygon into 2 or more polylines, in which case we no longer have an inside/outside region. During calculation of inverse transformation using iterative algorithm like NTv2, a null point can cause failure to compute the inverse even if the result would be perfectly valid because the point could be temporarily in the invalid region before the offset applied by iterative algorithm brings it back to the valid region. So in summary, the difficulty is not to tell that a point is null, but the consequences on geometries, rasters, inverse transformations, etc. The implications can be difficult enough to create a temptation for developers to fill null with their own values. For example in the case of points in geometry, there is a risk that developers will just ignore the null point and let the video card do a linear interpolation from previous valid point to next valid point.

Continuity on grid boundaries

From a developer perspective I think that continuity on boundaries of adjoining, overlapping, or nested grids is highly desirable. However I do not know if it would put too much constraints on producer side.

It has been pointed that bicubic interpolations can require the values of nodes from both grids near the edge. Actually even bilinear interpolations could require both grids, depending how the grids are defined: if defined like the left side in the figure in the "Continuity at grid border" section, even bilinear interpolations require both grids. The definition illustrated on the right side aims to avoid the need for two grids in the bilinear interpolation case, at the cost of an overlap (the green nodes). So supporting bicubic interpolations with a single grid may be seen as increasing an overlap that we already accepted for supporting the bilinear interpolation case.

What is a base grid

I would said a grid which is not a child (nested grid) of any other grid.

If there are overlapping base grids how should the correct value be defined

There is two scenarios. The simplest one is transforming points in isolation (i.e. not part of geometry, raster, etc.). In this case I think that PROJ convention is to iterate over the grids in a preference order defined by the producer, and to select the first grid containing the point? That convention seems fine to me.

The second and more tricky scenario is when transforming a geometry or raster. In this case I suggest to take the first (in producer-defined preference order) grid that contains fully the geometry or the raster. If no such grid is found (e.g. if the geometry always overlap at least 2 grids), then this is an open question…

Nested grids

About whether nested grids add or replace parent values, I was thinking to replacements for performance reasons. The addition strategy is of course doable too, but it seems to me that it has performance cost (2 bilinear interpolations during forward transformations, 2 iterations during inverse transformations). I would suggest to support at least the replacement strategy, and keep the addition strategy in option if desired.

Hello Kevin. I would said that it is the intention that GGXF allows the grid shifts to be applied directly in projected space. At least I see no reason to restrict to geographic space from an implementation point of view. Actually the space could be any CRS representable by ISO 19162 (or other CRS encoding) including parametric CRS, engineering CRS, etc. I think that the choice of restrictions to put on the coordinate spaces is more a conceptual problem than an implementation problem.

About the rules in Coordinate system axes section, the proposal is to prevent mixture of projected and geographic space between interpolation CRS and target CRS. But mixtures would be allowed between source CRS and interpolation CRS. Or to be more accurate the complete picture would be as below (arrows represent a coordinate conversion or transformation; I put a prime on second "interpolation CRS" because I'm not sure if it should be considered a separated CRS or not).

(Source CRS) ⟶ (Interpolation CRS) ⟶ (Interpolation CRS)′ ⟶ (target CRS)

The transformation using values interpolated on the grid is the arrow between (Interpolation CRS) and (Interpolation CRS)′. I propose that the axis restrictions (same direction, same units) apply between those two CRS. I'm trying to capture the idea that transformation is computed by taking source coordinates, applying a small translation on them and obtaining target coordinates that are only slightly different than the source coordinates.

All other arrows can be any ISO 19111 operation and GGXF would put no restriction on them. Consequently there would be no restriction (other than those implied by ISO 19111) between (source CRS) and (Interpolation CRS) axes. The intent is to allow the target CRS to be in a different space than the source CRS, but nevertheless keep that space change as a separated step.

Note that (Interpolation CRS) and (Interpolation CRS)′ are mostly conceptual CRS for discussion purpose; I'm not sure there is a need to formally express them, except maybe in transformations like IGNF (France) where source and target CRS are geographic, and interpolation CRS is geocentric.

Returning null, NaN or throwing an exception is easy when transforming points without worrying about their connections with other points. But a difficulty happens when the point is part of a geometry.

@desruisseaux I'm not sure this is any harder - we are basically saying that the whole geometry cannot be transformed to/from in a particular coordinate system (if that is what a specific GGXF is used for). Depending on the usage of the GGXF the producer could provide (or GGXF metadata could provide) a definition of an area of validity that users (or software) can apply trim to if this is what they want to do. I am not envisaging random null values scattered through the grid - more likely at the edges beyond the jurisdiction of the producer agency. And as you say the producer can choose to provide an interpolation if they feel that meets their users' needs better.

Also in our use cases for GGXF we include hydroid models. I'm not sure it is reasonable to define these inland? However I think the main thing is that we need to check our use cases and/or users and check that there isn't a desire or requirement for supporting null values.

@kevinmkelly I think that affine transformation in the description is the mathematical transformation of coordinates from the coordinate space (geographic or projected) to the GGXF grid coordinates (0,0), (0,1) etc. Would it be clearer to describe this in terms of ordinates (eg x,y) that could be longitude/latitude, easting/northing, or anything else for that matter.

True, if the recommendation is to declare the whole geometry as invalid when at least one point is null, then the problem become as easy as the simple point case. It is also true that declaring inland values in a hydroid model would be problematic. In the case of an hydroid model however, if I'm not misunderstanding, the values have no impact on the horizontal coordinates of geometries and rasters. In such case, at least for mapping on a flat horizontal map, null values are not problematic. This case is a little bit different than the "datum shift" case where values change the coordinates in the same dimensions (the horizontal ones) than the dimensions used for locating nodes on the grid.

At the very least my previous suggestion to avoid null values must become "it depends"…

Source and target CRS should have axes in the same order and same direction

In EPSG most geoid models are exposed as a transformation between a Geographic 3D CRS and a Vertical CRS, so violating this constraint. But such transformation can also be seen as between a Geographic 3D CRS and a compound CRS made of the Geographic 2D CRS of the source and the vertical CRS. But such compound CRS may not have a EPSG code, wheras its 2 components have one. Or there might be transformations between 2 VerticalCRS that use a Geographic 2D CRS as interpolation.

I think that PROJ convention is to iterate over the grids in a preference order defined by the producer

Actually what was ultimately implemented is to use the hiearchy of grids defined either explicitly by subgrids pointing to their parent by name (prefered way) or by the relationship of bounding boxes (in the absence of explicit relationships between grids). So when transforming a point, PROJ starts at the level of base grids. When it finds a base grid containing the point, it uses it, unless this base grids has children, in which case it repeats the process recursively. Normally, the base grids shouldn't intersect, as well as child grids at the same level of a same grid, but if they do, then which one is selected is implementation defined.

Regarding nodata, some existing grids do have nodata (thinking to Austria' geoid model for example). I believe we must support that, even if it is annoying for some applications.

In EPSG most geoid models are exposed as a transformation between a Geographic 3D CRS and a Vertical CRS, so violating this constraint.

Yes, the geographic 3D CRS is the source CRS, but I think we can describe that as having an implicit Geographic 2D interpolation CRS. In this case that interpolation CRS stay the same before and after the transformation.

It is true that some interpolation CRS does not have an explicit EPSG code, but above discussion about interpolation CRS is more an attempt to define the concepts. I'm not sure that the interpolation CRS needs to be explicitly identified, since it can sometime be implicit.

@ccrook said

Continuity ... Do we specify that it must be for a valid GGXF grid, or do we allow the possibility of discontinuity and provide metadata to inform software if the grid is continuous.

I suggest that there are issues here that we need to separate. Whilst from the perspective of fast access to any part of a grid and unambiguity in interpolation it is highly desirable that a grid is continuous and fully populated, I don't think we can prescribe this because there are going to be cases where it is impossible to achieve. I suggest we need three things: i) The format should allow for both continuity and discontinuity in the data, for extrapolated values to fill a 'grid block' and for a no-value flag; ii) there should be guidance for file producers on the advantages arising from them creating files that are 'recommended' - continuous, fully populated, etc; iii) there should be guidance for implementers on how best to treat files where the content has discontinuities, no values, etc.

Whether (ii) and (iii) are in the same document or because the audiences differ are separate documents, and whether this is(these are) in an informative annex in the format document or a separate guide, needs to be discussed.

Personally I favour a single format document including informative guidance so that everything is together, but it may depend upon how voluminous the guidance is to be.

desruisseaux said:

I'm not sure that the interpolation CRS needs to be explicitly identified, since it can sometime be implicit.

Yes it can be implicit, but when this is the case (e.g. in geoid models where it is implicit that the interpolation CRS is the horizontal 2D subset of the source geographic 3D CRS) this is unambiguous only when required by the operation method documentation. EPSG has decided that all new methods using grid interpolation should have the interpolation CRS explicitly given (as an operation parameter). Given that including method descriptions in the header is probably not going to happen because of the complexity of doing so in a machine-readable manner, to remove any potential for ambiguity, I believe that GGXF should make it a requirement in the file header for explicit declaration of the interpolation CRS, through either of the means described in 19111 clause 7.2: (a) a full description, (b) reference to a full description in a register of geodetic parameters, or (c) both (a) and (b),

with (b) recommended (as in 19111) but (a) is necessary to cover cases where the CRS is not in a registry. If (c) is to be permitted in a GGXF file, the 19111 conflict resolution statement "If there is a conflict between the two, the object full description should prevail over the reference to a register" should be included.

We need to agree on how the full description in (a) should be required in a GGXF file. WKT2 is an obvious candidate, but JSON following some yet-to-be-developed CRS-JSON schema is a possible alternative.

@RogerLott

Continuity ... Do we specify that it must be for a valid GGXF grid, or do we allow the possibility of discontinuity and provide metadata to inform software if the grid is continuous.

I suggest that there are issues here that we need to separate. Whilst from the perspective of fast access to any part of a grid and unambiguity in interpolation it is highly desirable that a grid is continuous and fully populated, I don't think we can prescribe this because there are going to be cases where it is impossible to achieve. I suggest we need three things:

Note: The discussion on continuity is provided in much greater detail in https://github.com/opengeospatial/CRS-Gridded-Geodetic-data-eXchange-Format/issues/7 (which I promised at the last meeting).

I'm not sure how you are thinking that continuity or missing values have any impact on speed of access to the grid. They relate to what is at the grid node, not how you retrieve it. There is a potential for ambiguity if the standard does not define how to handle missing values.

i) The format should allow for both continuity and discontinuity in the data, for extrapolated values to fill a 'grid block' and for a no-value flag;

I am very uncomfortable with the extrapolating across missing values (if that is what you are saying). The whole point of a missing value is to specify that the spatial function represented by the grid is not defined at some locations. If it is appropriate to extrapolate across this then it should be the producers responsibility to do this and provide an unambiguous full grid. I can't see any reason why GGXF would allow missing values, and then recommend or support extrapolation across them.

That doesn't mean a user (or software) can't extrapolate - it just means that if they do so they are not compliant with GGXF.

If GGXF does support missing values then it could be useful for the grid metadata to specify whether or not the grid is fully populated. GGXF software could use this to determine whether or not it can handle the grid (if it does not support missing values).

If GGXF does support missing values then it could be useful for the grid metadata to specify whether or not the grid is fully populated

I'd say that if a grid has at least one grid node set at the nodata value, it should declare a metadata item giving the nodata value. If a grid has no values at nodata, then this metadata item should normally be absent.

@RogerLott It all looks very comprehensive.

For the use of offsets in various dimensions, there is some prior and common practice in meteorology where the various combinations of offsets have been named as Arakawa-A to Arakawa-E.

Your current text seems to allow offsets that are not half-grid lengths. If this is intended, an example would be useful. If only half-grid lengths are envisaged, it should be stated much more clearly.

Offsets do have the edge cases of points north and south of the north and south poles respectively, but the maths all works out fine. (Lat/long points not on the earth are just complex numbers.)

An important point is that vector quantities at a pole with multiple grid points can be handled correctly too, using the maths, rather than a programming fudge.

HTH, Chris

Allow me an observation regarding the recommendation to deal with null values and extrapolation: The ‘domain of validity’ logic above in my view excludes any extrapolation as it would fall outside the ‘area of validity’. Similarly, any undefined (‘’null’, void’) node would automatically create a zone of ‘invalidity’ around it that ends only at the next ‘valid’ point. In other words, interpolation may only be allowed where the closest four neighbouring nodes are valid and a respective method is stated. (The best solution would be of course to re-calculate the value using the original spatial function.)

Much useful background in this issue so referenced on wiki page. Closing as not a current issue for specification.

opengeospatial / CRS-Gridded-Geodetic-data-eXchange-Format