3.1.1.  Discussion on keeping the Coordinate Systems simple

An ISO 19111 coordinate system is only a set of axes at no particular location, together with geometry formulas that are implied by the coordinate system class (not to be confused with the coordinate system unions, which are not covered in this section). For example, the use of a CartesianCS implies that angles, distances, surfaces, volumes, etc. are computed by formulas derived from Euclidean geometry. The use of a CylindricalCS or SphericalCS implies different sets of formulas for computing geometric properties. But all CoordinateSystem classes defined by ISO 19111 imply relatively simple and well-known formulas. This simplicity is necessary for performing tasks such as computing intersections or solving integral equations for computing surfaces. The most complex coordinate system class defined by ISO 19111 is the EllipsoidalCS. In that latter case, solving the integral equations requires series expansions (Clause 2.19), as done in map projections or geodesic (Clause 2.9) calculations. Extensions to ISO 19111 such as MinkowskiCS should not be much more complex.

Figure 2 — ISO 19111 Coordinate Systems

The simplicity mandated by the approach defined in the above paragraph may seem restrictive, but is an unavoidable step applied even with complex topographies. Complexity is usually managed by splitting a surface into small cells. Inside each cell, the “simple” mathematics of a coordinate system is assumed to be a good approximation. For example, the figure below zooms in on a single cell of a complex surface. If a software package computes the shortest path (the geodesic) between two points inside that cell by applying linear interpolations between (x₁, y₁, z₁) and (x₂, y₂, z₂) coordinate tuples, then that software is implicitly using a Cartesian coordinate system. The cells size (grid resolution) is chosen in such a way that approximating the real surface by a flat surface inside the cell area is considered an acceptable error.

Figure 3 — Surface with local Cartesian CS in each cell

In the above example, the shortest path inside a cell can be computed using the global (Clause 2.10) Cartesian coordinate system for the whole grid (gray axes), or a Cartesian coordinate system which is local (Clause 2.13) to the cell (red axes). In the example, both coordinate systems produce the same results, but the local CS is simpler because it can be reduced to a two-dimensional coordinate system (i.e., the perpendicular axis can be omitted) if all points lie on the surface. Coordinate transformations between those two systems are not discussed in this section, as these transformations require Coordinate Reference Systems (CRS) to provide context for the Coordinate Systems (CS).

The local coordinate system does not need to be of the same class as the global one. For example, the local CS could be an EllipsoidalCS if that is a better approximation of the real geometry inside the area of that cell, as illustrated in the figure below. In this case, doing geodesic calculations in the local CS produces different results than the cases illustrated in Figure 3.

Figure 4 — Surface with local ellipsoidal CS in each cell

In the above examples, calculations must be limited to a single cell. If the feature of interest spans more than one cell, then calculations must be performed separately in each cell and the results combined in some way. ISO 19111 does not specify how multi-cell calculations should be done.

An application can choose to use the same CoordinateSystem (CS) for all cells, even if axis directions vary slightly between different cells. Using a common CS means that all cells use the same geometric formulas with axes having the same units of measurement and the same orientations relative to the cell. The precise orientations of the axes are encapsulated in the Reference Frame, which is not part of a Coordinate System (it is part of a Coordinate Reference System). Using a uniform Coordinate System makes combining the results of adjacent cells easier, for example, by stating that all cells use a CartesianCS with measurements in meters. Therefore, volumes or areas computed in adjacent cells can be summed without unit conversions. This is true even for EllipsoidalCS using ellipsoids of different sizes and shapes. In the ISO 19111 model, the shape of the ellipsoid is not a property of the Coordinate System. Rather, it is a property of the Geodetic Reference Frame which, together with the EllipsoidalCS, forms a Geodetic Coordinate Reference System.

3.1.2.  Discussion on CRS domain of validity versus CS simplicity

An ISO 19111 Coordinate Reference System (CRS) is associated with exactly one Coordinate System (CS), except Compound CRS, which is omitted for the following discussion. A CRS is also associated, through subclasses of the datum class, with information about the CRS origin and “absolute” axis orientations (e.g., relative to stars). Consequently each cell in a grid such as Figure 3 or Figure 4 has its own local CRS because each cell has a different reference frame for different origin and “absolute” axis directions.

Coordinates can be transformed from the global CRS to the local CRS of a particular cell before computing geometric properties in a given cell. Note that in this scenario, the complexity of the irregular surface appears in the process of transforming coordinates from the global CRS to a local CRS. This complexity does not appear in any coordinate system before or after the transformation. In ISO 19111 modeling, the complexity is in CoordinateOperation, not in CoordinateSystem.

Figure 5 — ISO 19111 Operations and Coordinate Reference Systems

Given the restriction that each ISO 19111 CRS (excluding CompoundCRS) shall be associated with exactly one CS, and given that a CS should be mathematically simple, the CRS in Figure 3 and Figure 4 can be based on one of the following choices.

1. The global three-dimensional Cartesian coordinate system, in which case the complex surface plays no role. That surface is a feature inside the CS domain, but does not define the CS itself and is not involved in any calculation between pairs of (x, y, z) coordinate tuples.

2. A global two-dimensional coordinate system resulting from a curve-fitting of the complex surface. The result can be the surface of a plane, cylinder, sphere, ellipsoid, or something else defined outside ISO 19111. That approximation can be very good in some cases (e.g., Earth surface approximated by an ellipsoid). This option is nevertheless not discussed further in this Engineering Report.

3. A local two-dimensional coordinate system on the complex surface, i.e., the CS of one of the cells. That coordinate system may be Cartesian or ellipsoidal (among others) as illustrated in the two above examples, depending on the desired compromise between accuracy and mathematical simplicity.

Each alternative has pros and cons. Option #1 allows the referencing of any point anywhere in the domain covered by the global CRS. However, it does not allow accurate geodesic calculations (e.g., shortest path on the surface) between pairs of points. Conversely, option #3 is limited to a small area, but allows much more accurate geodesic calculations in that area.

Note also that complex numerical models, such as those based on fluid mechanic (meteorological and oceanographical models) or based on the equations of Einstein’s General Relativity (Annex E.2), do not change the fundamental approach. Those numerical models still need to split a large space into a grid (Clause 2.11) of cells small enough to allow the use of “simple” mathematics inside each cell with acceptable errors. The grids may be more sophisticated than grids of equilateral cells, but they are still grids with cells delimited by analytic functions.

The above discussion implies that if geometric properties or geodesics paths are required with high accuracy, then the associated calculations may require a CRS restricted to a local area. However, the area can be made larger by a choice of coordinate system class, or by a choice of kind of acceptable errors. This is partially illustrated for one-dimensional coordinate systems by the figure below. However, the discussion should be thought of as if the coordinate systems were two- or three-dimensional. The left side shows the case when Cartesian coordinate systems are wanted. No CartesianCS can be a good fit for the whole mountain, so the calculations must be restricted to smaller areas. The calculations can only be done between pairs of points on the same side of the mountain. If a calculation needs to be done between two points that are on each side of the mountain, it needs to be split into two parts, computed independently on each side, then combined somehow. In the ISO 19111 model, those two computations are done in two different Coordinate Reference Systems. However, if the coordinate system is based on a quadratic curve instead of straight lines, as shown in the right side of the figure below, then it becomes possible to do computation on the whole mountain with a single CRS. Therfore, computation between points on different sides of the mountain is no longer a problem.

Figure 6 — More complex coordinate system covering larger area

In a two-dimensional coordinate system, the use of EllipsoidalCS instead of CartesianCS allows covering a larger surface with a continuous function in a way similar to Figure 6. In the case of the Earth, using EllipsoidalCS is often an approximation good enough for covering the whole planet with a single CRS, such as WGS84. However, it does not change the general principle that EllipsoidalCS is still an approximation. On some irregular bodies, that approximation may not be good enough to support the use of a single EllipsoidalCS for the whole body with the desired accuracy.

3.1.3.  Coordinate System versus Coordinate Transformation roles in accuracy

A coordinate system based on a cubic curve would be more accurate or would cover a yet larger area than the examples shown in Figure 6, but at the cost of more mathematical complexity. Because the equations for intersections, surfaces, etc. can quickly become unsolvable, that mathematical complexity cannot be increased indefinitely. For example, ellipsoidal equations are only of degree 2, and already very difficult to solve. To cover a larger area, an alternative to increasing mathematical complexity is to keep a simpler coordinate system, such as CartesianCS, but accept some errors. For example, it is possible to favor accuracy in surface calculations at the expanse of accuracy in angles, or conversely. Controlling those errors is the art of map projections.

The overall accuracy does not depend only on the accuracy of formulas in the chosen coordinate system. Overall accuracy also depends on the accuracy of the transformation used for expressing points in that coordinate system. For example, Figure 7 below illustrates the results of a numerical experiment done at a point located at 45°N on Earth. The following test was repeated at distances varying from 1 meter to 1000 kilometers around the above-cited point.

• A circle is drawn around the center point (0°E 45°N) with a radius equal to the distance to test. The circle is determined by computing the locations of 360 points (each point is separated by 1°) using geodesic formulas on an ellipsoid published by Karney (2013).

• Distances and azimuths (horizontal angles from north) are computed for paths from the center point to the above-cited 360 points using the following methods.

  • Using SphericalCS formulas (distances are the arc lengths):

    • On an authalic sphere (sphere of the same surface as the ellipsoid) using the geographic coordinates as-is; and

    • On the same authalic sphere, but with geographic coordinates converted from the geodetic latitudes on the flattened sphere (Earth) to the geocentric latitudes on the authalic sphere. Those latitudes differ by about 21 kilometers. This operation is identified as a “datum shift” in Figure 7.

  • Using CartesianCS formulas (distances computed by the Pythagorean formula):

    • After projection to “Mercator (variant B)” with standard parallel at 45°; and

    • After projection to “Lambert Conic Conformal (2SP)” with the two standard parallels at 45°.

• The above calculations are compared with the expected distances and azimuths, which are the values used for drawing the circle in the first step.

The Y axis below shows the relative error, as a percentage of the distance shown on the X axis. For example, an error of 0.25% for a distance of 10 kilometers is an error of 25 meters. The error bars correspond to one standard deviation.

Figure 7 — Relative errors of distances around 0°E 45°N computed with different methods

An observation is that naive application of spherical formulas on Earth, with WGS84 latitudes and longitudes used directly, is not very good compared to other approaches. The relative errors are constant at around 0.13% over the range of tested distances, and the standard deviation is almost as large. A more sophisticated calculation converting the latitude values before applying spherical formulas gives better results, with relative errors at around 0.06% and negligible standard deviations. However, and maybe counter-intuitively, the use of CartesianCS gives better results than SphericalCS over small distances. Two factors contribute to that accuracy:

• the use of a “standard parallel” projection parameter configures the map projection for optimal accuracy around the latitude used in this test, namely 45°N; and

• the Mercator and Lambert Conic Conformal map projection formulas take into account the flattened shape of Earth.

On distances small enough, the geometrical distortions caused by using CartesianCS are smaller than the errors caused by ignoring Earth flatness in spherical formulas or caused by using a global sphere radius which may not be optimal for the test area. (Note: a sphere with geocentric radius computed at 45°N was also tested, but the results are not shown because they were similar to the authalic sphere case). In this test, the use of “Mercator (variant B)” projection is more accurate than spherical formulas for distances up to 100 kilometers and the use of “Lambert Conic Conformal (2SP)” is still more accurate even at 1000 kilometers.

The results of using “Popular Visualization Pseudo Mercator” are not shown in this graph because the errors are around 41% and would not fit in the Y axis scale. Those errors are caused mostly by the scale factor not well suited to 45°N. Other factors contributing to the errors are the facts that this projection ignores Earth flatness and uses the equatorial radius instead of some globally better value such as the authalic radius (radius of the authalic sphere).

3.1.4.  Implications for MinkowskiCS

MinkowskiCS is a four-dimensional coordinate system developed for use in Special Relativity. It ignores the spacetime curvature caused by gravity according to General Relativity (Annex E.2). However, even in the context of General Relativity, the Minkowski coordinate system can be useful as an approximation good enough at short distances from the object of interest. Such an approximation is like the use of CartesianCS in map projections, which are flat approximations of a curved surface considered good enough inside a known region.

Because Cartesian and Minkowski coordinate systems are flat spacetime approximations, it is tempting to conclude that geodesic (Clause 2.9) calculations performed in those coordinate systems can only be less accurate than geodesic calculations performed in a coordinate system taking the curvatures in account, such as a spherical coordinate system. But as shown in Figure 7, this is not necessarily the case because the overall accuracy depends on two factors: the fitness of the coordinate system to the local geometry and the accuracy of the coordinate operation applied to get there (Clause 3.1.3).

The choice of a coordinate system impacts the complexity of a fundamental object of study in General Relativity. In the context of Einstein’s relativity, spacetime geometry is described by a metric tensor. The metric can be represented as a matrix and defines how distances, angles, and other geometry properties are computed. While those metrics are of particular importance in relativity, this concept can be used in classical geometry as well. The Euclidean metric used in CartesianCS is the simplest:

Let p be the coordinate tuple of a point in space or space-time, and p_i (or p_j) be one of the coordinates of that tuple. For a small displacement dp, the small distance ds is given by:

For the Euclidean metric, the development gives the familiar Pythagorean formula:

By comparison, the metric at the surface of a SphericalCS is as below, assuming a radius of 1, coordinates in (latitude, longitude) order and units in radians. The main thing to note is that this metric is not constant. The matrix varies at every position, because it depends on the latitude φ. Consequently, computing the geodesic distance s² on a non-short path requires resolving an integral equation (not done in this ER).

MinkowskiCS can be seen as CartesianCS generalized to the spatiotemporal domain. Its metric is shown below, using the (- + + +) sign convention. Coordinates are in (t, x, y, z) order with the time coordinate already multiplied by c, as required by the MinkowskiCS definition proposed in this ER:

In the same way as the metric on the surface of a SphericalCS may be considered of manageable complexity, there are some alternatives to MinkowskiCS which are also of manageable complexity. For example, assuming a system with spherical symmetry residing in a vacuum, then the Schwarzschild metric (Annex E.7) may provide a computationally feasible solution for ds² in Formula (2). This ER does not propose classes for modeling those alternatives yet, but it is expected that if MinkowskiCS is accepted, other classes could follow the path.

MinkowskiCS may not provide sufficient accuracy in some situations in the same way as CartesianCS on the Earth’s surface is not accurate when the area of interest becomes too large. The accuracy on Earth can be improved by “upgrading” CartesianCS to EllipsoidalCS. Likewise, the accuracy in a relativist context may be improved by “upgrading” MinkowskiCS to another relativist coordinate system using a more complex metric. However, those “upgrades” cannot be such that the geometry formulas become unsolvable. When the combination of a large region and high accuracy exceeds what can be solved, splitting the region in smaller cells as illustrated in Figure 3 become unavoidable.

This Engineering Report does not explore how to combine calculations done in adjacent cells. Doing so may require the participation of other OGC working groups such as Discrete Global Grid Systems (DGGS). However, those combinations would probably require coordinate transformations between cells. This is where the complexity of geoid surfaces or gravitational fields are involved. But this complexity is in the transformation between two CRSs, not in the CRSs themselves.

3.1.5.  Engineering CRS for spacecraft

The coordinate reference system of a spacecraft can be described by an EngineeringCRS. This type of CRS was discussed in OGC 22-038r2 and an example is provided with the DART.xml file. ISO 19111 restricts the set of CS classes that can be associated to an EngineeringCRS. The allowed CS classes are affine, Cartesian, cylindrical, linear, ordinal, polar, and spherical implying that MinkowskiCS is also accepted because it is a subtype of CartesianCS. However, this is not necessarily the case for other metrics such as Schwarzschild (Annex E.7). Relaxing this restriction on CS type is proposed in Clause 3.2.3.3.

3.2.1.  Additions to OGC Topic 2/ISO 19111

Most classes listed below are new and do not require modifications to the ISO 19111 Standard, as the classes can be specified as an extension in a separate document. An exception is CelestialBody which may need to be included in ISO 19111 if an association from Datum to CelestialBody is added.

3.2.2.  Additions to OGC WKT 2 / ISO 19162

This section proposes new keywords for the Well Known Text (WKT) 2 format. Those additions require changes in the definitions provided in Backus-Naur form (BNF). The following table shows a summary of new WKT keywords. Those new keywords should be added (except Minkowski) to the Well-known text representation of coordinate reference systems Standard (OGC 18-010r11) §6.6 — Reserved keywords, with “TBD” replaced by the actual section numbers where the keywords are defined.

Table 5 — Summary of new WKT keywords

[ { <wkt separator> <alias> } ] ...

<alias>         ::= <alias keyword> <left delimiter> <alias text> <right delimiter>
<alias keyword> ::= ALIAS
<alias text>    ::= <quoted Latin text>

<celestialBody>          ::= <celestial body keyword> <left delimiter>
                             <celestial body name> <right delimiter>
<celestial body keyword> ::= CELESTIALBODY
<celestial body name>    ::= <quoted Latin text>

<geodetic reference frame keyword> <left delimiter> <datum name> <wkt separator>
[ <wkt separator> <celestialBody> ] <ellipsoid>  [ <wkt separator> <datum anchor> ]
[ { <wkt separator> <identifier> } ]…  <right delimiter>
[ { <wkt separator> <prime meridian> } ]

3.2.3.  Changes to ISO 19111

This section describes proposed changes to ISO 19111 abstract model. If accepted, those changes would require the publication of an ISO 19111 amendment and/or an OGC corrigendum.

3.2.3.2.  Remove ObjectUsage for generalization

In the ISO 19111 model, most classes inherit a set of properties from the IdentifiedObject parent class. Those properties are name, alias, identifier, and remarks. But there is also one additional property, namely domain, which is inherited only by the Datum, CRS, and CoordinateOperation classes. ISO 19111:2019 does this restriction by defining those properties in an intermediate class named ObjectUsage. ISO 19111:2007 did this restriction by repeating properties in the three above-cited classes, in a way similar to the UML diagram provided in section Clause 3.1.

The following UML diagram shows the ISO 19111:2019 model. The ObjectUsage contains only a domain property and is inserted between the IdentifiedObject parent class and the Datum, CRS, and CoordinateOperation subclasses. This insertion was done for preserving the same restriction as in ISO 19111:2007, where only Datum, CRS, and CoordinateOperation could have a scope and domain of validity.

Figure 12 — Current ISO 19111:2019 model with object usage

Removing the restriction saying that only Datum, CRS, and CoordinateOperation can have a scope or a domain of validity is proposed. There is nothing inherently wrong in specifying a scope on most kinds of CRS-related objects. For example:

• A MinkowskiCS definition may want to specify “for use with special relativity” as the scope.

• An Ellipsoid definition may want to specify a limited geographic domain of validity if the celestial body has an irregular shape with different ellipsoidal approximations for different geographic areas.

• A PrimeMeridian definition may want to specify a limited temporal extent if it is determined by an ephemeral feature, for example of the surface of a gas planet.

• A CoordinateSystemAxis with direction such as “South along 60°E meridian” may want to specify ”for use at North pole” as the scope.

The restriction can be removed by deleting the ObjectUsage class (not to be confused with the ObjectDomain class, which is kept unchanged), and moving the domain attribute up to the IdentifiedObject parent class.

3.3.1.  Chain of operations

ISO 19111 coordinate operations can be chained together in a ConcatenatedOperation object, which contains:

• the overall source CRS of the chain;

• the overall target CRS of the chain;

• each step as a reference to a single operation (in GML documents, it can be a xlink:href attribute); and

• other attributes common to all operations (identifier, name, scope, domain, and accuracy).

Each step other than pass-through can be qualified as a Conversion, Transformation, or PointMotionOperation. The difference between conversion and transformation was discussed in the Testbed-18 Engineering Report OGC 22-038r2 §5.6, and illustrated in the Testbed-18 3D+ Data Space Object Engineering Report OGC 23-011r1 §8 — Importance of datum. The key takeaway is that transformations are the steps where stochastic errors are introduced. Thus, by checking the type of each step, the ISO 19111 model enables users to perform some analysis of the (possibly many) sources of errors. For example in Figure 14, there are stochastic errors in the step represented by a red box. Different chains of operations between the same (source, target) pair of CRS may use different transformation steps with different sources of errors.

A concatenated operation can have as many steps as desired, where each step may use a different method. The main constraint is that the target CRS of a step must be the source CRS of the next step, except for inverse operations where the above-cited CRS can be swapped. Specifically, if a concatenated operation is defined as below:

• Overall source: CRS A

• Overall target: CRS C

• Step 1: from CRS B to CRS A

• Step 2: from CRS B to CRS C

then the software should understand that step 1 needs to be inverted, i.e., it needs to go from CRS A to CRS B. This is illustrated by Figure 15. Sometimes, the software can derive the inverse operation simply by inverting the sign of parameter values. Sometimes, the inverse operation requires completely different formulas (e.g., map projections). How to derive the inverse operation shall be specified by the operation method documentation. This flexibility avoids the need to define every OperationMethod twice: once for the forward direction and once for the reverse direction. Users may define an explicit OperationMethod for the reverse direction if that is wished, but it is not required. The figure below is the case where it is not desired.

Figure 15 — Example of inverse operation in a concatenated operation

3.3.2.  Operation encoding

The Well-Known Text (WKT) 1 format cannot encode coordinate operations, except to some extent via the PARAM_MT keyword of the OpenGIS Coordinate Transformation Service Implementation Specification OGC 01-009, but the latter was not a well-known feature. The fact that WKT 1 was practically restricted to Coordinate Reference System encoding may have contributed to the perception that ISO 19111 provides no coordinate operation framework. OGC 18-010r11 (a.k.a. WKT 2) can encode some coordinate operations, but does not cover the full set of ISO 19111 features (e.g., it excludes PassThroughOperation) and has no mechanism such as GML xlink:href for avoiding repetition (e.g., each targetCRS is repeated in the sourceCRS of the next step). The Geographic Markup Language (GML) is the most capable encoding for coordinate operations in current OGC Standards, except for the new ISO 19111 classes introduced after 2007. Unfortunately, few open source software applications provide good GML support for coordinate operations. At the time this ER was written, Apache Spatial Information System (SIS) is the only open source project to our knowledge with this capability, but this project is not as widespread, for instance, as GDAL.

OGC 22-038r2 demonstrated the use of GML for describing a chain of coordinate operations from an imaginary observatory on Earth to an imaginary spacecraft. However, due to technical constraints of the prototype at that time, all steps and their CRSs had to be described in a single GML document, resulting in a XML file about 700 lines long. In this ER, the main components are described in separate GML files, with relationships specified by xlink:href attributes in GML. For example, a concatenated operation can be described as below:

<gml:ConcatenatedOperation gml:id="DimorphosToDART">
  <gml:identifier codeSpace="TB19-D001">DimorphosToDART</gml:identifier>
  <gml:name>Dimorphos to DART</gml:name>
  <gml:scope>Testbed-19 demonstration.</gml:scope>
  <gml:operationVersion>1</gml:operationVersion>
  <gml:sourceCRS xlink:href="BarycenterToDimorphos.xml#DimorphosCRS"/>
  <gml:targetCRS xlink:href="BarycenterToDART.xml#DART"/>
  <!--
    The first step below should actually be "Dimorphos to barycenter", but the software
    should detect, by inspection of the source and target CRS, that it needs to use the
    inverse of that operation.
  -->
  <gml:coordOperation xlink:href="BarycenterToDimorphos.xml"/>
  <gml:coordOperation xlink:href="BarycenterToDART.xml"/>
</gml:ConcatenatedOperation>

Figure 16 — Concatenated operation example in GML

The use of those files by a prototype is discussed in Clause 3.4.

3.3.3.  Finding relationships between CRSs

A measurement may be expressed in the local (proper) reference system (e.g., GeodeticCRS) or in that of an observing entity (e.g., EngineeringCRS). The value of a measurement may be expressed in an unlimited number of reference systems. ISO 19111 provides a framework to support coordinate transformations between arbitrary coordinate reference systems. The process can proceed as follows.

1. For each pair of Coordinate Reference Systems of interest, it describes the relationship between those two reference systems with a CoordinateOperation. An unlimited number of CoordinateOperation instances can be created, including many instances for the same pair of reference systems if desired.

2. If chaining some of the above-defined operations is useful, describe those chains with ConcatenatedOperation instances. It is not necessary (neither feasible) to enumerate all possible chains. The purpose of explicit ConcatenatedOperation instances is to state that, among all possible sequences of operations going from CRS A to CRS C, the path A → B → C (for example) is preferred.

3. List all above definitions in a database. For example, a relational database, or a directory containing a GML file such as Figure 16 for each coordinate operation. The EPSG geodetic parameter dataset is an example of a relational database defining a set of coordinate operations (in addition to the CRS definitions).

4. Implement the ISO 19111 RegisterOperations class, in particular the following method:

findCoordinateOperations(source : CRS, target : CRS) : Set<CoordinateOperation>

Figure 17 — Search of Coordinate operations in a register

The above method first searches in the database, or in the directory of GML files, for coordinate operations defined for the given source and target reference systems. If at least one CoordinateOperation is found, then all matching operations are returned. Users can inspect the metadata listed in Clause 3.3 (e.g., domain of validity and positional accuracy) for deciding which operation is best suited to the need. Users can also base decisions on an analysis of operations steps, for example as illustrated in Figure 14.

If no predefined operation was found for the given pair of reference systems, then the implementation can try to infer a sequence of operations. For example, if the user requested an operation from A to C but the database contains only an operation from A to B and another operation from B to C, then the implementation can infer the A → B → C sequence. However, this approach is not always deterministic as the number of possible combinations may grow exponentially with the number of steps, and the implementation may not be able to afford the cost of exploring all possible paths. Providing explicit ConcatenatedOperation instances in the database for preferred operation chains is desirable.

3.3.4.  Executing coordinate transformation

Once a CoordinateOperation has been chosen, coordinates can be changed from being referenced to the sourceCRS to being referenced to the targetCRS, and/or (in the case of dynamic CRS) from being referenced to one coordinate epoch to being referenced to a second coordinate epoch, by executing the following method.

transform(points : CoordinateSet) : CoordinateSet

Figure 18 — Coordinate operation execution

A CoordinateSet contains coordinate metadata (CRS and data epoch), together with an ordered list of coordinate tuples as DirectPosition objects. The transform method preserves the tuples ordering.

3.3.5.  Standard method proposals

The ISO 19111 framework for coordinate operations can be applied in space with no change. New operation methods need to be defined, but this is done using the existing model. For example, seven operation methods were proposed in OGC 22-038r2 annex A:

• four methods for expressing GeoPose’s Basic-YPR (yaw, pitch, and roll), Basic-Quaternion, Regular-Series, and Irregular-Series functionalities in terms of ISO 19111 model; and

• three methods for implementing a “Voyager to Observatory” scenario defined in that ER.

This ER focuses on only one operation method but defines that method more extensively for making it more suitable to standardization. The operation method name is “Translation by trajectory” and the proposed code in OGC namespace is http://www.opengis.net/def/method/OGC/1/200 (this URL would also be the place for information on where to publish the definition). The definition is specified as a GML encoding in Annex B.2. This first operation method is very simple. It:

• works on 4-dimensional (x, y, z, t) coordinate tuples, but non-relativist;

• assumes that axis directions relative to distant stars do not change; and

• transformations consist of translations provided by a Moving Feature file.

Note that the simplicity of the above characteristics is not due to ISO 19111 limitations. There is no constraint in the number of parameters that an operation method can define. Those parameters can be Booleans, integers, list of integers, measures, list of measures, character strings, filenames, URLs, citations, or geometries. Parameters can also be defined in groups and repeated. The OGC 22-038r2 ER used parameter grouping to demonstrate how a piecewise function can be defined with ISO 19111 framework. In this Testbed-19 ER, the operation method defines only one parameter:

• http://www.opengis.net/def/parameter/OGC/1/200 — “Trajectory file”

The value of that parameter is a filename or URL to a Moving Feature file in a format defined by OGC: JSON, XML, CSV, or netCDF (for this Testbed prototype, only CSV is supported). The Moving Feature file defines the trajectory of the target CRS relative to the source CRS. The target CRS may be, for example, an EngineeringCRS attached to a spacecraft and the source CRS may be an InertialCRS attached to a planet. The CRS declared in the Moving Feature file for trajectory data is taken as the interpolation CRS. In the simplest case, the interpolation CRS is the same as the source CRS. The algorithm to apply is as follows.

1. For each (x, y , z , t) source coordinate tuples, interpolate the position in the trajectory at time t. The result is taken as a translation vector in units of the interpolation CRS.

2. If the interpolation CRS is not the same as the source CRS, then:

   • get the derivative (Jacobian matrix) of the transform from source CRS to interpolation CRS at the location given by source coordinates;

   • inverse the Jacobian matrix; and

   • multiply the inverted matrix by the translation vector obtained in step 1. The result is the translation vector converted from interpolation CRS to source CRS.

3. Subtract the translation vector from the source (x, y , z) coordinates and pass the value t unchanged.

The algorithm for the reverse operation is below. The condition in step 2 should not mention the source CRS, but in the case where the difference between source and target CRS is only a translation with no scale or rotation, the derivative of the transformation between those two CRSs is an identity matrix everywhere and step 2 can be skipped.

1. For each (x, y, z, t) target coordinate tuples, interpolate the position in the trajectory at time t. The result is taken as a translation vector in units of the interpolation CRS.

2. If the interpolation CRS is not the same as the source or target CRS, then:

   • get the derivative (Jacobian matrix) of the transform from target CRS to interpolation CRS at the location given by target coordinates;

   • inverse the Jacobian matrix; and

   • multiply the inverted matrix by the translation vector obtained in step 1. The result is the translation vector converted from interpolation CRS to target CRS.

3. Add the translation vector to the target (x, y , z) coordinates and pass the value t unchanged.

The use of this operation method is demonstrated by the prototype in Clause 3.4.

3.4.1.  Implementation-specific interpretations of OGC Standards

There are two issues that are not formally addressed by OGC standards (decisions left to implementers), but which must be resolved to produce an executable demonstration (Annex B.3). The first issue is in which directory are the CSV files listed in Table 6? This ambiguity is not specific to Testbed-19. For example, the EPSG database also references datum shift files (NADCON, NTv2, etc.) in the same way as this ER references the Moving Feature files. Implementations such as Apache SIS or PROJ expect those datum shift files to pre-exist in dedicated directories. But this approach may not be applicable for transformations using user-supplied trajectory files. For the Testbed-19 prototype (Annex B.3), it has been decided that the “Trajectory file” parameter value (e.g., as in Figure B.4) is interpreted as a path relative to the location of the GML file.

The second issue is how to resolve the CRS identifiers in Moving Feature files. The various encodings for Moving Feature (MF) do not include full CRS definitions, but only use identifiers (e.g., Figure B.5). But because of the space context, those identifiers are not EPSG identifiers. In Testbed-19, the identifiers are rather like urn:ogc:def:crs:JPL::120065803 (note: the use of “JPL” in this ER does not mean an endorsement of the Jet Propulsion Laboratory. It is only an illustration of what it may look like if JPL was publishing those definitions). Since there is no CRS registry for the “JPL” authority, where should users look for this information? The prototype developed and documented in this ER takes advantage of the fact that full CRS definitions are provided in GML documents, as part of the ISO 19111 CoordinateOperation definition. For example, the BarycenterToDimorphos.csv MF file (Figure B.5) is referenced by the BarycenterToDimorphos.xml GML file (Figure B.4). The latter itself contains xlink:href attributes to full CRS definitions like Figure B.2, which contains the following line:

<gml:identifier codeSpace="JPL:HORIZONS">urn:ogc:def:crs:JPL::120065803</gml:identifier>

Figure 19 — Fragment of Didymos CRS definition

The Apache SIS library uses this context for tracing urn:ogc:def:crs:JPL::120065803 back to its full CRS definition. Note that there is no need to explicitly create a project-specific registry in some file. Instead, a “virtual registry” is created in memory only and populated automatically by the library at GML parsing time.

3.5.1.  CRS registry

There is currently no Celestial Reference System (CRS) such as BCRS (Barycentric Celestial Reference System) in the OGC registry. OGC maintains an http://www.opengis.net/def/crs/IAU/ namespace for astronomical bodies, but this namespace currently contains only Coordinate Reference Systems for planetology based on the work of the IAU’s Cartographic Coordinates & Rotational Elements working group. Allocating identifiers for CRSs from other sources is not straightforward, as the OGC Planetary DWG bases the identifiers on NAIF identifiers for uniqueness and semantic reasons.

However, the absence of Celestial Reference Systems in the OGC registry is not necessarily a blocking issue. Clause 3.4 demonstrates how celestial reference systems can nevertheless be used, with CRS definitions in GML files referenced by chains of coordinate operations defined in other GML files. Formats that do not include full CRS definitions, such as Moving Features (MF) datasets, can nevertheless reference those definitions by their identifiers when the MF files are read in the larger context of reading the GML files. Therefor, having inertial CRS definitions in a authoritative registry is a convenience but is not mandatory, as the same CRSs can be defined outside any authoritative registry and still be found by the software.

3.5.2.  General relativity

General Relativity may be handled with a spatiotemporal geometry modeled like terrain. It may be possible to leverage the techniques developed for the 3D geoid for 4D spacetime. However, at the time this report was written, the geoid data formats and tools are disparate. The first attempt to provide a unified format as an international standard is the Gridded Geodetic data eXchange Format (GGXF). Testbed-19 participants did not experiment with the possibility of using GGXF because the GGXF Standard is new and its support in software is still in development.

4.4.1.  Separate Model(s) and Encoding(s)

GeoPose 1.0 has both conceptual and logical models and a corresponding JSON encoding. This packaging came about primarily because the market forces pushing GeoPose standardization expected a JSON realization. The quickest path to a usable standard was to specifically target JSON, and to combine models and encoding in a single document.

Clarity of design and reusability of building blocks within the design is facilitated by separation of models and encoding. Derivation of one or more logical models from a conceptual model or at least a body of well-defined concepts also encourages compatibility and interoperability of the components. NeoPose prototypes a more nuanced approach to a logical implementation of GeoPose concepts and interfaces to supporting elements such as coordinate systems and coordinate reference systems defined in a way that facilitates using existing libraries for implementation of frame transformations.

4.4.2.  Fully Defaulted Optional Properties

The extensions to GeoPose in this ER are described and implemented by adding new properties to the GeoPose conceptual model. These new properties potentially add to the size and complexity of GeoPose for the most common applications where a small amount of information is required. This size burden can be avoided by requiring that the corresponding standard specifies well defined default values for each of these new properties. If a computation or other use of a GeoPose requires a property that is not explicit in a serialization, that unique default value is used. Thus a value for each of these is always available, though not necessarily explicitly in a serialized form.

4.4.3.  Uncertainty

For this experiment, a simplified approach to uncertainty was used.

Quantitative parameters were given an optional precision property in the same unit of measure as the parameter itself. This precision property is a symmetrical bound on the range of values within a 90% confidence limit.

The precision of the result of a coordinate transformation is also important. In some cases it is possible to quantify how imprecision in input parameters results in imprecision in the transformation results. This imprecision in the transformed coordinates is a combination of the transformation itself and imprecision that may be introduced by an approximate method, e.g., linear regression or approximation by series expansion (Clause 2.19). NeoPose does not consider this an important type of error.

Qualitative parameters were given a fractional belief property in the interval [0,1]. A value of 0 indicates no belief that the qualitative category is correct. A value of 1 indicates that the value is certainly correct and intermediate values indicate intermediate belief that the value is correct.

4.4.4.  Time

Since the standardization of GeoPose 1.0 in 2022, the importance of time has become more obvious with experience since the standardization of GeoPose 1.0 in 2022. Especially for poses involving moving entities, time is a key element to synchronize poses derived from independent sensor systems. The concept is the simple one of a 1-dimensional time coordinate axis and the location of points — instants — along that axis. A temporal coordinate reference system is needed in order to establish the time coordinate axis, in line with OGC proposals for the modeling of time in OWL time. A time coordinate axis is a one-dimension coordinate system counting time from a starting instant, or epoch.

In TB-19, the time coordinate axes are either UNIX Time or GPS Time. The default value is 0. The default precision is 0.000000000001 second (one picosecond), representing 0.3 mm of light travel time in a vacuum. This value was chosen to be compatible with the precision of terrestrial surveying measurements.

4.4.5.  Identity

In most applications, a GeoPose is a property of an entity in the physical world. Associating the entity and GeoPose can be accomplished by modeling the association as a property of an object representing the entity, as a property of the GeoPose, or as an independent bidirectional association. In the case where the association is a property of the GeoPose, there are additional considerations.

First, in GeoPoses representing sensed and detected entities, there can be uncertainty regarding both the existence of the entity and its globally unique identity. For example, a traffic camera may detect a car but initially not detect a license plate. Assigning the detected car a temporary identity with a short “time-to-live” solves the need for a locally unique identifier. This can be used for purposes such as tracking persistence of semantic type through time, as well as temporary hiding of part or all of a detected entity behind a closer entity. It is also useful in resolving changing classifications of a single physical entity through time.

The TB-19 solution was to add an optional unique id and an optional explicit time-to-expire property.

4.4.6.  JSON Serialization

Extensions should support a JSON serialization that is as backwards compatible as possible with OGC GeoPose 1.0. Backwards compatibility still allows definition of a logical model that differs from GeoPose 1.0 as long as the JSON encoding still conforms to the JSON-Schema definitions.

4.4.7.  Support for Operations

Implementation of the GeoPose Standard with reference frames outside the default 3D WGS-84 geodetic coordinates must make it easy to use existing libraries (e.g., PROJ) or to recognize the frames from parameters (e.g., ISO 19111) or IDs (e.g., EPSG numeric identifiers).

An updated GeoPose logical model should better support implementations of frame transformations defined by existing external standards. A simple serializable structure for Outer and Inner frames, as well as the specific transformation between them, would make it possible for implementers to encapsulate the details. This would better support the use of existing libraries as well as new implementations in a single structure. This framework for modular implementations would not expose the complexity of the multiple forks of some commonly used libraries.

4.5.1.  Architecture

The overall idea is that there is a first transformation that converts the coordinates of an event $e_O=\left(\begin{array}{c}t\\ x\\ y\\ z\end{array}\right)$ in the outer frame to corresponding coordinates in an inner frame: $e_I=\left(\begin{array}{c}{t}^{\prime }\\ x^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)=T_{OI}\left(e_O\right)$

To be useful, the Outer frame should be attached directly, or via a sequence of transforms, in a linked GeoPose structure to an externally-defined coordinate reference system. This gives the event $e_I$ in the Inner frame context in that external system of time and position. This makes the pose “Geo” if the external reference is tied to the Earth.

The second transform defines a ray cast from the origin of the inner frame. The second transform can take many forms but it is also a transformation of an event $e_R=\left(\begin{array}{c}{t}_R\\ x_R\\ y_R\\ z_R\end{array}\right)=T_R(e^{\prime },s)$, where s is a real number parameter, e.g., distance along the ray in the length coordinates of the Inner frame. A ray may be defined in many ways, some more or less suitable depending on the application context. For example, the ray can be defined as a sequence of rotations, or as a quaternion, or as direction cosines, or as a point on the ray not at the origin of the inner frame CS, or as a rotation matrix, or as another transformation that defines a ray (or geodesic, to be more general) – such as a Lorentz transformation.

These transformations are free to assume default values for unspecified input coordinate values and to ignore, i.e., to project transformed coordinates. For example, time may be ignored, or a 1- or 2-dimensional space may be used.

To summarize, the core structures in the NeoPose areas follows.

• The PoseID which gives a pose a permanent unique ID.

• Pose links which are lists of predecessor and successor poses’ poseIDs.

• A frame transform, including specification of inner and outer frames, as well as a transformation with any needed parameters.

• A further transformation within the inner frame, often a rotation.

The second transformation represents the orientation element of the two-part pose paradigm: position plus orientation. When any pose — not only GeoPose — is represented by data, the orientation is defined by the parameters of a coordinate transformation that is often rotation or rotation-like, such as the Lorentz transformation. These parameters are in addition to the coordinates themselves, which are along a ray (or more generally, a geodesic) rooted at the origin of the inner frame coordinate system.

Both transformations take a time coordinate value and three spatial coordinate values (i.e., an event) to a transformed event. The first transform redefines an event in an external CRS with a description of that event in the inner CRS. It is the same event. The second transform provides the ray or geodesic, which is a path within the inner CRS. This is a set of events, parameterized by an independent variable such as distance that generates the ray or geodesic.

Many transformations come from families of transformations, which are identifiable by sets of controlling parameters. Examples include the WGS-84 coordinates in the first transform of each of the GeoPose 1.0 Basic Forms, the quaternion of the Basic-Quaternion form, and the yaw, pitch, and roll angles in the Basic-YPR form.

Figure 22 — The NeoPose Core

4.5.1.1.  Idiomatic and Templated Transforms

An event is a point in space-time. This point has three spatial and one temporal coordinate. Every NeoPose transform maps an event, that is, three space coordinates and one time coordinate, to another event. In many cases events may be projected, with one or more of the coordinates ignored in the transform.

The specification of transforms in a NeoPose can be of two types:

• idiomatic (identified by parameters); and

• templated (identified by an externally standardized scheme).

This division is the NeoPose strategy for dealing with the wide ranges of frame transform specifications across organizations and disciplines. An analogy to the issue of recognition and interpretation is recognition of a cat in an ordinary photo vs in a stylized symbol such as the toothy smile of the Cheshire Cat in Stevenson’s Alice in Wonderland.

Idiomatic transforms: The key issue in implementing frame transforms is the existence of multiple idioms for specifying the controlling parameters. Some are derived from computer graphics, some from mathematics, some from geodesy and cartography, and some from navigation. Most of these idioms are broadly used only within a limited range of use cases. It can be awkward and confusing to force a “natural” expression of a transform using a paradigm from another field. The NeoPose approach is to accept that there are useful idiomatic expressions of transforms tied to specific disciplines or use cases.

Templated transforms: In other cases, it is normal to use an organizing scheme (template), such as the concepts of ISO 19111, OGC WKTs, by reference to EPSG identifiers, SEDRIS SRM, or NASA SPICE kernels. NeoPose represents such templated transform expressions with key-value pairs. The keys are arbitrary but unique, and the values are expressions according to one of the standardized organizing schemes. No attempt is made within the NeoPose representation to parse or interpret these structures. Instead, a form of dynamic or late binding is used. The transformations within the used scheme must itself be able to parse the expressions and use them to implement the described transformation.

There is a separation between transforms that are fully defined by the OGC GeoPose standard and transforms that follow some other system of definition and operational application. There is a mechanism to use foreign or external specifications of reference frames and transformations between them. There is also a mechanism to bring the transformations into the OGC standard by defining associated properties as a new idiomatic transform.

Figure 23 — NeoPose Abstract Transform

4.5.1.3.  Concrete Forms

The concrete idiomatic transforms implemented for TB-19 are shown in the following figure. These include support for the GeoPose 1.0 Basic-YPR and Basic-Quaternion forms as well as a newly proposed “Local” form pair similar to the Basic forms except that the outer frame is a Cartesian 3D coordinate system with lengths in meters and time in picoseconds.

Figure 24 — NeoPose Concrete Transforms

The “GeodeticPosition” object in “WGS84-LTPENU Transform” in the above figure is the “position” object from the GeoPose 1.0 Standard Basic forms. The “YPRRotation” and “QuaternionRotation” objects are the corresponding “angles” and “quaternion” objects. These are examples of the selection of specific transforms from families using transform-specific selection parameters. The following figure has examples of the serialization of the selection parameters.

The following figure shows some examples of serialized forms showing the parameters selecting specific transforms from several idioms from GeoPose 1.0 as used in TB-19. Again, the serialized parameters are not input to the implemented event transforms but rather are the method of selection of specific member transforms sharing an idiom.

Figure 25 — Serialized NeoPose Idiomatic Transforms

Although templated transforms could be implemented with serialization conforming to GeoPose 1.0, this was not done for TB-19. The following is an implementation of the WGS84 to LTP-ENU transform previously shown as an idiomatic transform. Every idiomatic transform can be serialized as a templated transform. Identifying a serialized templated transform as an idiomatic transform is a difficult task and not generally useful. Using the templated form is a reasonable means of documenting idiomatic forms.

Figure 26 — NeoPose Templated Transforms

4.5.1.4.  NeoPose Support Datatypes

The following are four datatypes that are used in the NeoPose implementation.

Figure 27 — NeoPose Datatypes

The property NeoContext is a convenience element supporting the use of long text strings, such as OGC WKTs defining coordinate reference systems or transformation. The object holds key-value pairs where a short key (local to and defined by an application, system, or protocol) can be used to refer to a long, complex string. For TB-19, a table of 14 entries were implemented as a WKT JSON “virtual registry”.

{
  "position": {
    "lat": 47.7,
    "lon": -122.3,
    "h": 11.5
  },
  "quaternion": {
    "x": 0.2303923729814588,
    "y": -0.03749706362805052,
    "z": 0.15620176329098345,
    "w": -0.9597470155390087
  }
}

4.5.2.  Implementing Selected Suggestions

The extensions developed in TB-19 are intended to adopt some of the suggestions from TB-18 and to demonstrate that these changes also provide support for GeoPose in non-Euclidean spaces. The name “NeoPose” is used for the purpose of distinguishing the extended GeoPose from the version specified in the OGC 1.0 Standard.

In possible future work, some or all of the NeoPose extensions may also be standardized.

4.7.1.  Class hierarchy

One way to compare NeoPose and ISO 19111 models may be that NeoPose sub-typing of Transform (with the idiomatic and templated transform subclasses — Figure 23) addresses encoding brevity aspects. ISO 19111 sub-typing of CoordinateOperation (with conversion and transformation subclasses among others — Figure 5) addresses more conceptual aspects about the nature of coordinate operations (e.g., whether the parameter values are “by definition” or stochastic). The representation of a transform in ISO 19111 model provides rich metadata but without trying to be compact. Conversely, the representation of the same transform in NeoPose model is more compact but loses some metadata.

4.7.2.  Transform parameters encoding

NeoPose uses parameter names directly in a JSON serialization. For example in Figure 25, the “Cartesian-Translation Transform” uses “dx,” “dy,” “dz,” and “dt” parameter labels, while the “Cartesian-YPRRotation Transform” uses “yaw,” “pitch,” and “roll” parameter labels. By comparison, ISO 19111 parameters are also named, but in an indirect way. Instead of:

{
    "lat": 47.7,
    "lon": -122.3,
}

Figure 29 — Snippet of NeoPose idiomatic transform parameters

The ISO equivalent (for illustration only, as there is no OGC Standard yet for JSON encoding of ISO 19111):

{
    "parameterValue": {
        "parameter": "lat",
        "value": 47.7
    },
    "parameterValue": {
        "parameter": "lon",
        "value": -122.3
    }
}

Figure 30 — Snippet of ISO 19111 equivalence (non-standard)

The description, unit of measurement, range of values, etc. of “lat” and “lon” parameters are defined in a separated class. That separated class can be seen as describing the template. This topic was covered in more detail, together with examples of GeoPose transform definitions, in OGC 22-038r2 annex A.3.1.

4.7.3.  Transforming positions or events

NeoPose provides the following method (Figure 23):

TemplatedTransform.transform(outerCRS : CRS, innerCRS : CRS, Parameters, Event) : Event

Figure 31 — NeoPose transform method

By contrast, ISO 19111 separates the process into the following two steps (Clause 3.3.3).

RegisterOperations.findCoordinateOperations(source : CRS, target : CRS) : Set<CoordinateOperation>

CoordinateOperation.transform(points : CoordinateSet) : CoordinateSet

Figure 32 — ISO 19111 transform method

The first step (findCoordinateOperations) is typically much more costly than the second step (transform), because it may require a search for all coordinate operations registered in a database for the given pair of source and target CRS. By contrast, the second step may be as simple as subtraction. Because the first step typically needs to be executed only once before to transform thousands of points, there is a big performance advantage in separating those two steps. It allows executing the costly part outside a loop, then to loop only on the cheapest part. Another reason for this separation is that findCoordinateOperations may return many transforms, so the user needs to choose one before to moving to the next step.

However, the performance argument is to be nuanced. The NeoPose transform method expects a Parameters argument, which is absent from the ISO 19111 methods because those parameters are to be found by findCoordinateOperations in the database. Providing the parameters may reduce the needs for NeoPose to perform database queries. An inconvenience is that it puts the burden on users for finding the parameter values. An advantage is that it is easier for users to change those parameter values if desired.

4.7.4.  Contextual registry

NeoPose defines a NeoContext class used as a contextual registry (Clause 4.5.1.4). This registry associates CRS definitions to identifiers using key-value pairs in a file. This file must contain all CRSs referenced by NeoPose transforms (Table 7).

By contrast, the prototype based on the ISO 19111 model takes a different approach (Clause 3.4.1 and Annex B.3). It also has a contextual registry but built automatically in memory. There is no file listing explicitly all key-value pairs. The prototype takes advantage of the fact that ISO 19111 definitions of CoordinateOperation can contain full CRS definitions, even if indirectly through xlink:href attributes. The Apache SIS library used in the prototype automatically builds an in-memory virtual registry populated with all CRS definitions encountered during GML parsing and uses that registry for resolving CRS identifiers in Moving Feature files referenced by the CoordinateOperation. However, this strategy is an implementation choice. The OGC/ISO standards make no recommendation.

4.7.5.  NeoPose features not covered by ISO 19111

Both NeoPose and ISO 19111 can associate identifiers to reference systems. However, NeoPose has the concept of temporary identifiers with an optional explicit time-to-expire property. There is no standard equivalence in ISO 19115 (the metadata standard used by ISO 19111 for identifiers). However, implementation-specific equivalences are possible by creating subclasses of MD_Identifier.

NeoPose transforms not only positions, but also orientations. ISO 19111 is only about positions, at least in the text of that Standard. However, all classes and methods operate on coordinates that are only identified by their position in coordinate tuples. It is technically possible to use 6-dimensional coordinate tuples (or 7-dimensional if including time) with position in the first 3 coordinates and orientation in the last 3 coordinates, then to use those tuples in coordinate operations. Whether this would be stretching too much the spirit of ISO 19111 has not been debated.

4.7.6.  ISO 19111 features not covered by NeoPose

ISO 19111 provides more metadata than NeoPose. ISO 19111 uses ISO 19115 (metadata) data structures for expressing citation, identifiers, spatial, vertical, and temporal extents, and ISO 19157 (data quality) for expressing uncertainty. Each of those Standards go far in expressing, for example, how the uncertainty was estimated. Then ISO 19111 itself has a rich set of classes for expressing the semantic of coordinate operations (e.g., conversion versus transformation) and how the operations are chained together (e.g., Figure 14). NeoPose does not yet go as far, except for the orientation aspect which is not addressed by ISO 19111. However, ISO 19111’s richness makes it more difficult to encode.

5.2.1.  Barycentric Celestial Reference System (BCRS)

The BCRS has its origin at the Solar System’s center of mass, and its axes are aligned with the International Celestial Reference Frame (ICRF). The ICRF is fixed relative to far away extragalactic objects, defined by measuring quasar positions.

The IAU defined it as a system of space-time coordinates for the solar system in the year 2000. It is considered the system appropriate for the basic ephemerides of solar system objects and astrometric reference data.

The use of the BCRS and of the Barycentric Coordinate Time (TCB, from the French temps coordonnée barycentrique) provides a single consistent frame of reference minimizing the effects of General Relativity for observations within our Solar System.

The reference plane of the BCRS is the celestial equator, whereas the primary direction (reference longitude) is the March (Vernal) equinox. BCRS coordinates can be expressed in 3D Cartesian form (x, y, and z) in units such as Astronomical Units (au) or kilometers (km), as used in the JPL Horizons input data and the MF-JSON files generated by these experiments. The BCRS coordinates could alternatively be expressed in spherical coordinates, using declination (δ or lat), right ascension (α or lon), and range coordinates.

Figure 34 illustrates the axes names, reference plane, and primary direction for the Geocentric Celestial Reference System (GCRS), which are the same as for the BCRS. For the BCRS, the Earth reference center would simply be replaced by the Solar System Barycenter.

Figure 34 — ICRF Axes for GCRS (source: Wikipedia)

5.2.2.  Barycentric Coordinate Time

The Barycentric Coordinate Time (TCB, from the French temps coordonnée barycentrique) is the time coordinate complementing BCRS.

TCB is equivalent to the proper time for a clock at rest in the coordinate frame, moving together with the Solar System center of mass, but not influenced by the gravity of solar system objects such as the Sun.

TCB replaces the older Barycentric Dynamical Time (TDB, from the French temps dynamique barycentrique).

6.1.1.  The Hillyfields Bubble Data Capture

Ultimately, the purpose of many spatial and spatiotemporal representations is to act as a replica of a bounded part of the physical world. A space-time bubble (“bubble”) is one type of replica. A bubble replica is a bounded unit of spatiotemporal information. Within the bubble’s bounds, there are representations of all physical objects and materials relevant to the intended application of the replica. Paper maps are an example of a bubble. A paper map (or the map in your car’s navigation system) has well-defined bounds, a period of temporal validity, mapping of semantics to symbology, and anything important according to its modeling criteria is represented. In addition, the cartographic representation must be interpretable and useful to humans. To generalize a bit, a replica (e.g., a bubble) must:

• represent objects bounded in space and time, i.e., the bubble surface;

• have a semantic representation;

• be complete under specific conditions, i.e., meet a “closed-world” condition (this means that nothing exists but the contents and that the contents completely fill the time and space bounded by the bubble surface); and

• provide access to the contained semantic elements anchored in space and time in a way that supports the retrieval and manipulation requirements of a user experience.

The bubble concept is very useful in consideration of GeoPose objects and especially in composite collections of related GeoPose objects.

One way to define a bubble is with a GeoPose and a set of bounds. To illustrate use of OGC GeoPose 1.0 in complex configurations and to explore extensions of GeoPose to support a completely general pose concept that works in any type of space or space-time, TB-19 used a bubble representing a real location, centered on the headquarters of the UK OS. Using a single well-defined replica that contains a variety of semantic elements whose locations, orientations, semantic categories, and appearance were observed in several timespans. This bubble can serve as a Testbed for developing, testing, and comparing software and systems in interactive applications and experiences, especially those involving complex pose relationships between objects in multiple reference frames.

Figure 40 — Hillyfields Bubble, UK Ordnance Survey

The Hillyfields Bubble is a multipurpose dataset captured in support of the development of Augmented Reality/Metaverse architectures and software infrastructure prototypes. Data capture was successfully carried out at the UK OS on April 25th and 26th, 2023. The capture included data from visual and navigation sensors on:

• a Land vehicle [RTK GPS, compass, gyro, accelerometer, 2 video streams, LiDAR];

• two people [GPS, gyro, compass, accelerometer, audio, 5312×2988 video];

• two fixed video cameras recording with Web Video Map Tracks (VMT) annotation; snf

• a drone following car [GPS, compass, gyro,1920×1080 video], 7 runs 25th, 5 runs 26th, 10 Hz pose resolution.

There were eight runs on the 25th and three runs on the 26th.

6.1.2.  Simulated Satellite in Low Earth Orbit

Simulated coordinates were calculated for each second from a point overhead Longyearbyen, Norway to a point over Bamako, Mali.

JSON track or KML track with one point per second.

6.1.3.  Extension to Minkowski Space-time

Two tasks were defined in an area centered on the car park in front of the headquarters building of the UK Ordnance Survey in Southampton, UK:

• Pointing task: Aiming an imaging sensor in low earth orbit at a moving vehicle on the ground; and

• Minkowski space visualization task: Visually observing a collection of entities in approximately uniform relative motion at velocities approaching the speed of light.

These tasks were designed to exercise new features in the extended GeoPose 1.0 (NeoPose proposal) model. Data from instrumented moving vehicles, persons, and a fixed background environment were collected on April, 2023. This data capture in a constrained volume of space and time was named the “Hillyfields Bubble”.

6.2.1.  NeoPose Representation of Objects within the Hillyfields Bubble

Each frame of each imaging sensor attached to each moving object was processed to identify people, cars, cameras, and drones. The lists of observed objects were represented as Local NeoPoses. In some cases, it was possible to identify an external imaging sensor and to establish the identity of an observed object.

6.2.2.  Minkowski Space-time Visualization Task

In this task, each of the camera sensors may be the observer and the background environment as well as each of the other moving sensors are in corresponding observed frames. In Galilean relativity, the geometry of all objects in observed frames is viewed as though the speed of light were infinite. In other words, the geometries of the observed objects are independent of their motion relative to the observer. In Minkowski flat space-time, the view of observed objects must be determined by casting rays back from the projected pixels of the sensor. Although one might consider applying a Lorentz transformation to the entire geometry of an observed object in uniform relative motion, this “coordinate” view will not give a correct visual or sensor result for an extended object. What is necessary is, in effect, to project the projected image of the sensor image plane backwards in time to find what objects are intersected. If there are k objects in different inertial frames, then k versions of each traced ray must be considered individually, and the results combined to show the actual resulting light reaching the observer sensor pixel. The resulting computational cost is then k times that of classical raytracing in Galilean relativity. The actual ray casting and intersection computation is not changed:

1. For each pixel $p$ :

   1. Set the color $C_p$ to black and the number of contributing rays $N_p$ to 0:

      1. For each moving object $O_k$:

         1. cast a ray from the image plane through a camera pinhole;

         2. Lorentz transform the ray $R_{p_k}$ to object frame of $O_k$;

         3. find intersection, if any, with object k; and

         4. if there is an intersection, then using some (any) lighting model and object material, compute a color at the intersection point and update $C_p$ and $N_p$.

These steps must be repeated for every frame of an animation sequence. In the Hillyfields Bubble, this is a rate of 29.97 frames per second. For an image size of 1920 x 1080, and k=5 moving objects, this is 310 million rays per second. Raytracing takes several seconds per frame.

6.2.3.  Details for the Lorentz Transformation

Suppose the observed object is moving at a relative velocity $v^{ij}$ with respect to the observer:

Then compute

and

Then for object $i$

The actual transformation is

6.2.4.  The Pointing Task

Determination of the departure and arrival directions of a geodesic (Clause 2.9) between pairs of events in distinct inertial frames is a simple definition, but calculation of pointing angles, i.e., the orientation part of a pose, is in general not easy because an observer needs to have a target event in the frame of the source. There must be some source of predictions of the world line of the source in order to compute the geodesic from source to observer.

In TB-19, as the spacecraft orbit is visible above the horizon at the Hillyfields Bubble, the clock is set to 0 and a signal is sent to the Hillyfields Bubble with the location of the spacecraft in its frame using the Hillyfields LTP-ENU coordinate system, as well as a prediction of the future location each second for the next 500 seconds. A copy of the location message can be found here.

Figure 41 — Satellite Track from Longyearbyen, Svalbard to Bamako, Mali.

In Galilean or Minkowski space-time, geodesics are unique. They are also straight lines in Euclidean space. In general space-time, there may be many curved geodesics linking source and observer. That case is not considered in this ER.

Figure 42 — The Pointing Task in Space-time.

• Aim visual sensor on satellite in orbit at 750 000 km above the Hillyfields z = 0 plane at LiDAR unit on Street Drone:

  • about 140 seconds of visibility time;

  • approximate inertial reference frame for this time; and

  • compute pointing angles with respect to spacecraft frame.

The computations are the same as for raytracing, except that the rays are cast from the intended target back toward the satellite track. Intersections with the predicted satellite track give a sequence of rays that will intersect the desired object at some time during the satellite’s overpass.

7.4.1.  Vehicle Tracking

Moving objects were tracked from video using a process that was designed to take advantage of the oblique viewing angle by using lines of sight. This method is also robust to variations in the pan and tilt of the tripod mounting which inevitably occurred when the camera was removed to download files.

Figure 45 — Lines of Sight With Target Location Zones

The procedure can be summarized as follows:

1. identify suitable lines of sight from the video frames;

2. identify target locations corresponding to those lines;

3. identify video frames with a moving object that matches those target locations;

4. add a path segment to the corresponding WebVMT file; and

5. verify accuracy using previsualization in a web browser.

Figure 46 shows the observed track synchronized with video and a faint line indicating the camera’s line of sight.

Figure 46 — Vehicle Tracking Previsualization

7.4.2.  Synchronization Issues

The next challenge was to accurately synchronize data from different devices with a common timeline — in this case Coordinated Universal Time (UTC). Timing errors hinder the process of correctly associating observations of the same object at different times and locations between different devices. Hence, poor synchronization significantly reduces the value of data aggregation.

The smartphones used to record video from static tripods were linked to the cell network which enables access to accurate UTC time. However, smartphone cameras typically have a significant delay between a user request to start video recording and the first frame of that recording, and no accurate measurement of the video start time could be found. The latency depends on a number of factors including the camera hardware response time, resource pressure from other tasks, and the device’s power management strategy. As a result, the estimated start time was only accurate to within a few seconds.

By contrast, the StreetDrone vehicle was able to record timing data accurately to within 10ms and is synchronized with UTC. A process was devised to match the trajectory recorded aboard the vehicle with the trajectory observed by the smartphone cameras to synchronize their video timelines with UTC. Visual analysis is limited to the unit time of a single video frame which is around 33ms at 30 frames per second.

Data were recorded natively in WebVMT by the smartphones and exported to WebVMT from the StreetDrone’s inertial measurement unit (IMU) enabling tools developed in Testbed-17 to be reused for analysis and allowed results to be integrated with online maps for verification in a web browser.

A new WebVMT utility was developed to calculate a list of closest approach times to a static location by a moving object. This tool can calculate results to millisecond accuracy and includes a radius parameter that enables multiple approaches to be discriminated correctly. By matching the StreetDrone trajectory data against the target locations used for vehicle tracking by the roadside cameras, an average time offset was calculated for each video file captured by the static cameras. In addition, the standard deviation was calculated to quantify the quality of the synchronization process. Results were used to adjust the WebVMT start time parameter in order to accurately synchronize the path segments observed by the static camera with UTC.

A second utility was added to the existing suite that clips WebVMT content to match video files with revised start times and durations. WebVMT files synchronized with UTC were clipped and merged using this and an updated version of the Testbed-17 tools suite to create previsualizations that demonstrate the accuracy of these results — see Synchronized Previsualization.

Figure 47 — Synchronized Data Previsualization

7.5.1.  Synchronization Metrics

The geotagged video files recorded had fairly consistent time offsets from UTC for each of the two smartphone devices — one was typically 1.5 seconds ahead and the other was around 2 seconds behind UTC. Quality metrics were good for both cameras since most standard deviations are in the range between 50ms and 200ms which is indicative of the average error.

Table 9 — Nokia 5 Synchronization Metrics

Table 10 — Moto G30 Synchronization Metrics

The best results represent a statistical error of less than 100ms (3 frames) with a 95% confidence level. In terms of road vehicles, this equates to a distance error of less than 1.3m at 50km/h (30mph) or 3.1m at 110km/h (70mph).

The accuracy of these results is confirmed by the previsualization screenshots shown in Synchronized Previsualization.

7.5.2.  Wrong-Way Traffic Use Case

The scheme devised to track passing vehicles from video footage could be easily extended to detect wrong-way traffic.

Target location zones are sufficiently small to allow discrimination of individual lanes on a road. Using three or more zones per lane enables an automated detection routine to reliably determine in which direction a vehicle is traveling, and to discriminate between adjoining lanes which may have traffic flowing in opposite directions.

Figure 48 — Target Zones To Detect Wrong-Way Traffic

Figure 48 shows the target location zones and lines of sight in blue. Detection labels for the zones in lanes A and B are shown in orange, and the arrows show the correct direction of traffic flow.

The expected sequence of detections for traffic traveling in the correct direction is either A1→A2→A3 or B1→B2→B3→B4 in this case — since vehicles drive on the left in the UK. Wrong-way traffic would be detected by either A3→A2→A1 or B4→B3→B2→B1. This system could be easily adapted for use on a multi-lane road such as a motorway but could also be used to police a no-overtaking zone on a two-way road.

An unexpected result was discovered during analysis which demonstrates the strength of this traffic monitoring scheme. The StreetDrone performed a U-turn which was captured by one of the static cameras. Although the algorithm was not designed to identify this maneuver, the incident was detected since the detection zones were triggered in an erroneous order that correctly indicates the event and demonstrates the value of this simple approach.

Figure 49 — Unexpected U-Turn Detection

7.5.3.  Vehicle Speed Forensics Use Case

Another application for sight line analysis techniques was highlighted by an episode of the BBC Crash Detectives documentary about the Gwent Police Forensics Collision Investigation team — first broadcast on 9 October 2023.

In the first episode of the fourth series, the team investigated a hit and run incident in which a man died on a Welsh mountain road in the early hours of the morning. Although there were no witnesses or CCTV cameras, officers were made aware of video footage from a home security system located about 1km away on the other side of the valley. The cameras had recorded timestamped video which included car headlights traveling along the relevant road in the distance.

Initial police analysis of the night-time footage was unpromising, but then a daytime image was obtained from the same camera which they overlaid to show reference points in the scene. This approach allowed officers to estimate the vehicle’s speed at the time of the accident.

The problem seemed well-suited to the techniques developed during this Testbed, so the vehicle’s speed was estimated using WebVMT analysis in order to demonstrate how easily this could be applied to an operational police issue.

Figure 50 — Crash Detectives WebVMT Analysis

Two full-screen CCTV clips from the show were analyzed. The camera location was accurately identified to within 1m from freely available web maps to establish relevant sight lines and detection zones on the road.

Table 11 — Vehicle Speed Analysis

Results were calculated for two segments of the observed journey based on geospatial analysis of WebVMT data using the same sight line techniques previously described.

7.7.1.  Litter Monitoring Use Case

Away Team captured data for the litter monitoring use case in April 2023, which could be studied in Testbed-20 since there was insufficient time to address this in Testbed-19.

Litter locations were identified and exported to WebVMT which could be matched against video and LiDAR data recorded by the StreetDone. Analysis could be performed using WebVMT tools developed during Testbed-17 for tracking static objects from a moving vehicle to monitor litter locations.

Figure 51 — Identified Locations For Litter Monitoring Use Case

In addition, video was concurrently captured by StreetDrone’s front- and rear-facing cameras to enable multi-view data aggregation. Hence the risk of other road vehicles obstructing the survey vehicle’s view of a target location during a single pass can be mitigated in operational use.

Synchronizing data gathered by different devices and aggregating data from multiple passes of the same target location enables accretion rates to be calculated over time enabling more accurate monitoring and prediction of litter levels so maintenance teams can be deployed more efficiently.