3.11.  Abbreviated terms

CRS

Coordinate Reference System

CS

Coordinate System

CSV

Comma Separated Values

ECEF

Earth-Centered — Earth-Fixed

ECI

Earth-Centered Inertial

GGXF

Gridded Geodetic data eXchange Format

GML

Geographic Markup Language

WKT

Well-Known Text

1D

One-dimensional

3D

Three-dimensional

4D

Four-dimensional

5.1.  ISO 19111 overview

ISO 19111 supports mixing spatial, parametric, or temporal aspects in a single coordinate reference system. Such a coordinate reference system is called a compound CRS and the spatial, parametric, or temporal parts are components of the compound CRS.

A coordinate reference system (CRS) may be geodetic and apply to a national or regional scale. Such a coordinate reference system is called a geodetic CRS or projected CRS.

A coordinate reference system may also apply locally such as for a building or construction site. Such a coordinate reference system is called a derived CRS or engineering CRS, depending how they are defined.

Finally, a coordinate reference system may be referenced to a moving platform such as a car, a ship, an aircraft, or a spacecraft. Such a coordinate reference system is called an engineering CRS. The CRS of a spacecraft can be related to a second coordinate reference system which is referenced to the Earth through a transformation that includes a time element. This is discussed in more details in Clause 5.2.

The definition of a coordinate reference system does not change with time. In object-oriented programming, coordinate reference system objects shall be immutable. However, the defining parameters of a CRS may include a rate of change of the parameter. By this mechanism, and others, ISO 19111 supports the construction of CRS in which the coordinate values may change with time, for example due to tectonic plate motion. A reference frame with such time evolution is called a dynamic reference frame and should not be confused with a moving feature or dynamic feature. Both dynamic reference frame and moving feature describe the motion of points. But in the ISO 19111 context, a dynamic reference frame accounts for the deformations of the object associated to the frame, not for the displacement of that frame relative to another frame, such as the displacement of a moving feature relative to Earth. ISO 19111 currently restricts dynamic datums to geodetic and vertical reference frames. However, if the concept was generalized to spacecrafts, the motions of points in a dynamic datum would describe the deformation of the spacecraft itself. This idea will be explored in Clause 5.9. For now, spacecraft are considered as rigid bodies and ignore dynamic datums. In other words, only static (in ISO 19111 sense) datum are considered.

The following UML is taken from OGC Topic 2 and shows the main (not all) classes. The coordinate reference systems mentioned in the above paragraphs can be seen. The key elements to retain from this UML for the purposes of this ER are as follows.

• There is a base CRS class with many sub-classes specialized for different usages.

• All CRSs except CompoundCRSs are associated with a datum (or an ensemble of datum) and a coordinate system.

• DerivedCRSs are associated with a base CRS and a Conversion.

Figure 1 — UML of main ISO 19111 classes (source: OGC Topic 2)

5.2.  CRS referenced to a moving platform

For creating a coordinate reference system associated with a moving platform, an engineering datum is first defined. All ISO 19111 data contains an anchorDefinition field of type CharacterString. That field contains a description of the relationship used to anchor a coordinate system to the object. The anchor may be an identified physical point with the orientation defined relative to the object. The engineering datum itself is not time-dependent, but any transformations of the associated coordinates to an Earth-fixed or other coordinate reference system shall contain time-dependent parameters. An example is given in Clause 5.7.

In this ER, a CRS referenced to a hypothetical spacecraft as illustrated below is defined. This will be used as the source of a chain of coordinate transformations described step-by-step in subsequent sections. As a scenario, imagine that an astronomical body has been observed by the spacecraft. The body coordinates are relative to the spacecraft, using the CRS illustrated below. The goal is to transform those coordinates to a CRS relative to an observatory on Earth, using only existing OGC/ISO standards with as few extensions as possible.

The above CRS can be encoded in Well Known Text (WKT) format such as in the figure below. The text in the EngineeringCRS element can be anything and should be descriptive. The text in the Anchor element is mainly for human consumption. But the text in EngineeringDatum (i.e., the datum name) is significant. Implementations often use the datum name or identifier for deciding if two Datum elements are the same frame, because there is no other property in ISO 19111 datum that could be used. Note that despite a common practice, there is no need for a “D_” prefix or for replacing spaces by underscores.

ENGINEERINGCRS["A spacecraft-centred CRS for Voyager",
  ENGINEERINGDATUM["Voyager spacecraft frame",
    ANCHOR["Spacecraft centre of gravity"]],
  CS[Cartesian, 3],
    AXIS["Forward (x)",  forward],
    AXIS["Sideward (y)", starboard],
    AXIS["Downward (z)", down],
    LENGTHUNIT["metre", 1]]

Figure 3 — Engineering CRS for Voyager (WKT)

Following is an equivalent snippet in GML. This is a fragment of the VoyagerToObservatory.xml file in the GitHub repository. This fragment provides the same content as the above WKT, with the addition of scopes and identifiers (mandatory in GML), more explicit separation of axis names and abbreviations, and explicit name spaces for controlled vocabularies. The issue of GML verbosity compared to WKT is discussed in Annex B.

<gml:EngineeringCRS gml:id="VoyagerCRS">
  <gml:identifier codeSpace="TB18-D025">VoyagerCRS</gml:identifier>
  <gml:name>A spacecraft-centred CRS for Voyager</gml:name>
  <gml:scope>TestBed 18 demonstration.</gml:scope>
  <gml:cartesianCS>
    <gml:CartesianCS gml:id="VoyagerCS">
      <gml:identifier codeSpace="TB18-D025">VoyagerCS</gml:identifier>
      <gml:name>Spacecraft coordinate system</gml:name>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="VoyagerAxisForward">
          <gml:identifier codeSpace="TB18-D025">VoyagerAxisForward</gml:identifier>
          <gml:name>Forward</gml:name>
          <gml:axisAbbrev>x</gml:axisAbbrev>
          <gml:axisDirection codeSpace="ISO">forward</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:axis>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="VoyagerAxisStarboard">
          <gml:identifier codeSpace="TB18-D025">VoyagerAxisStarboard</gml:identifier>
          <gml:name>Sideward</gml:name>
          <gml:axisAbbrev>y</gml:axisAbbrev>
          <gml:axisDirection codeSpace="ISO">starboard</gml:axisDirection> <!-- Not well defined in this example. -->
        </gml:CoordinateSystemAxis>
      </gml:axis>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="VoyagerAxisDown">
          <gml:identifier codeSpace="TB18-D025">VoyagerAxisDown</gml:identifier>
          <gml:name>Downward</gml:name>
          <gml:axisAbbrev>z</gml:axisAbbrev>
          <gml:axisDirection codeSpace="OGC">down</gml:axisDirection> <!-- Not well defined in this example. -->
        </gml:CoordinateSystemAxis>
      </gml:axis>
    </gml:CartesianCS>
  </gml:cartesianCS>
  <gml:engineeringDatum>
    <gml:EngineeringDatum gml:id="CenterOfMass">
      <gml:identifier codeSpace="TB18-D025">VoyagerFrame</gml:identifier>
      <gml:name>Voyager spacecraft frame</gml:name>
      <gml:scope>TestBed 18 demonstration.</gml:scope>
      <gml:anchorDefinition>Spacecraft centre of gravity</gml:anchorDefinition>
    </gml:EngineeringDatum>
  </gml:engineeringDatum>
</gml:EngineeringCRS>

Figure 4 — Engineering CRS for Voyager (GML)

A CRS such as the one above does not contain relationship(s) with any other CRSs. The relationship between a spacecraft and the reference frame of a planet is encoded separately. The link uses ISO 19111 CoordinateOperation element and is described in Clause 5.6.

5.3.  Local tangent plane to ECEF

ISO 19111 supports the definition of new coordinate reference systems as the application of a coordinate conversion from another pre-existing coordinate reference system. Such a coordinate reference system is called a derived CRS, and the coordinate reference system from which it was derived is called the base CRS. A derived CRS inherits its datum or reference frame from its base CRS. This constraint is enforced by the use of a coordinate conversion instead of a coordinate transformation (see Clause 5.6) in the derived CRS definition. This implies that if a base CRS is a geodetic CRS referenced to Earth, then the derived CRS is still a geodetic CRS referenced to Earth. This construct cannot be used for transforming coordinates from a geodetic reference frame to the engineering reference frame of a spacecraft. However, derived CRS can be used for the definition of a local coordinate reference system tangent to the geoid or ellipsoid of a planet (among other possibilities).

The following WKT snippet describes the CRS illustrated above. When the METHOD element references an operation method defined by some authority like EPSG (which is the case in this example), the CRS definition contains enough information to allow map projection libraries to convert coordinates from Earth to the platform, and conversely. Operation methods are discussed in more details in Clause 5.3.1.

GEODETICCRS["Local Tangent Plane - East North Up (LTP-ENU)",
  BASEGEODCRS["WGS 84",
    DATUM["World Geodetic System 1984",
      ELLIPSOID["WGS 84", 6378137, 298.257223563]],
    ID["EPSG",4979]],
  DERIVINGCONVERSION["WGS84 to LTP-ENU",
    METHOD["Geographic/topocentric conversions", ID["EPSG",9837]],
    PARAMETER["Latitude of topocentric origin", 28.583333, ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",8834]],
    PARAMETER["Longitude of topocentric origin", -80.583056, ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",8835]],
    PARAMETER["Ellipsoidal height of topocentric origin", 100, LENGTHUNIT["metre", 1], ID["EPSG",8836]]],
  CS[Cartesian, 3],
  AXIS["Topocentric East (U)", east, ORDER[1], LENGTHUNIT["metre", 1]],
  AXIS["Topocentric North (V)", north, ORDER[2], LENGTHUNIT["metre", 1]],
  AXIS["Topocentric height (W)", up, ORDER[3], LENGTHUNIT["metre", 1]],
  USAGE[
    SCOPE["Voyager 2 launch"],
    AREA["Space Launch Complex 41"],
    BBOX[28, -81, 29, -80]]]

Figure 6 — Derived CRS for observatory platform (WKT)

5.3.1.  Operation method in derived CRS

A derived CRS contains (indirectly) an operation method (METHOD[…]). The operation method is a critical element of information which cannot be described fully in WKT or GML. “Method” is rather a controlled vocabulary for identifying a set of formulas published in some external document. A good example of such document is EPSG guidance note #7-2 for map projection formulas. Another example is the Pole rotation method defined by the OGC MetOcean working group. A third, more incomplete example is shown in Annex A.3. This example lists the operation methods derived from GeoPose or needed by VoyagerToObservatory.xml. That annex is much less complete than previous examples because it does not describe the formulas. That work would need to be done by an OGC Standard Working Group (SWG).

The formulas described in the above-cited documents need to be hard-coded using some programming language (C/C++, Java, etc.) in map projection libraries. Then the WKT and GML formats can reference those formulas by specifying the operation method name and identifier (”Geographic/topocentric conversions” and EPSG:9837 in above example). GML goes one step further than WKT by requiring a more detailed <gml:OperationMethod> element, which includes the list of parameters accepted by the operation method (only their descriptions, not the values, which are provided separately). The following UML is taken from OGC Topic 2 and shows the classes for describing an operation method. The key elements to retain from this UML for the purposes of this ER are as follows.

• OperationMethod has references to a formula and a set of parameters.

• Formula is either a human-readable text or a citation to some publication; there is no executable code.

• Description of parameters (GeneralOperationParameter) are separated from values (GeneralParameterValue).

Figure 7 — UML of operation method classes (source: OGC Topic 2)

The <gml:OperationMethod> element can be either embedded within the derived CRS definition or be provided elsewhere and referenced using xlink. The VoyagerToObservatory.xml file in the TB18 GitHub repository uses the embedded approach for implementing the full “Voyager to Observatory” scenario in a single standalone file. By contrast the examples in this ER use the xlink approach, which separates the operation method from the CRS definition. This separation is desirable because the operation method is typically defined by some authority such as EPSG and shared by an arbitrary number of derived CRS instances. For example, EPSG defines a few tens of operation methods. These methods are shared by thousands of projected CRS definitions in the EPSG database and countless WKT strings outside the EPSG database which rely on EPSG method definitions.

The first GML fragment (below) provides an operation method definition which is basically a simplified copy of a definition provided by the EPSG database. Note that formulas are specified only by the “See IOGP Guidance Note #7-2” text. A human must read that document for implementing the formulas.

<gml:OperationMethod gml:id="ECEF_to_LP">
  <gml:identifier codeSpace="IOGP">urn:ogc:def:method:EPSG::9837</gml:identifier>
  <gml:name codeSpace="EPSG">Geographic/topocentric conversions</gml:name>
  <gml:formula>See IOGP Guidance Note #7-2.</gml:formula>
  <gml:sourceDimensions>3</gml:sourceDimensions>
  <gml:targetDimensions>3</gml:targetDimensions>
  <gml:parameter>
    <gml:OperationParameter gml:id="epsg-parameter-8834">
      <gml:identifier codeSpace="IOGP">urn:ogc:def:parameter:EPSG::8834</gml:identifier>
      <gml:name>Latitude of topocentric origin</gml:name>
    </gml:OperationParameter>
  </gml:parameter>
  <gml:parameter>
    <gml:OperationParameter gml:id="epsg-parameter-8835">
      <gml:identifier codeSpace="IOGP">urn:ogc:def:parameter:EPSG::8835</gml:identifier>
      <gml:name>Longitude of topocentric origin</gml:name>
    </gml:OperationParameter>
  </gml:parameter>
  <gml:parameter>
    <gml:OperationParameter gml:id="epsg-parameter-8836">
      <gml:identifier codeSpace="IOGP">urn:ogc:def:parameter:EPSG::8836</gml:identifier>
      <gml:name>Ellipsoidal height of topocentric origin</gml:name>
    </gml:OperationParameter>
  </gml:parameter>
</gml:OperationMethod>

Figure 8 — Geographic/topocentric conversion method (GML)

The CRS definition can then reference the operation method definition using the xlink attribute. The following fragment assumes that the XML fragment shown in the above figure was in the same GML document. The ECEF_to_LP label is arbitrary as any name is okay because it is local to the document:

<gml:method xlink:href="#ECEF_to_LP"/>

Referencing definitions in external documents is also possible. For example, the above <gml:OperationMethod> fragment could be completely omitted from the XML document and the above link replaced by the following one:

<gml:method xlink:href="https://epsg.org/api/v1/CoordOperationMethod/9837/export?format=gml"/>

The XML fragment shown below uses such an external link for <gml:baseCRS> elements in order to keep the fragment smaller. However, the VoyagerToObservatory.xml file contains the full definition.

<gml:DerivedCRS gml:id="ObservatoryPlatformCRS">
  <gml:identifier codeSpace="TB18-D025">ObservatoryPlatformCRS</gml:identifier>
  <gml:name>Local Tangent Plane - East North Up (LTP-ENU)</gml:name>
  <gml:scope>TestBed 18 demonstration.</gml:scope>
  <gml:conversion>
    <gml:Conversion gml:id="WGS84_to_LP">
      <gml:identifier codeSpace="TB18-D025">WorldToPlatform</gml:identifier>
      <gml:name>WGS84 to LTP-ENU</gml:name>
      <gml:scope>TestBed 18 demonstration.</gml:scope>
      <gml:method xlink:href="#ECEF_to_LP"/>
      <!--
        This is where the parameter values for the conversion from ECEF to the local tangent plane are provided.
        Note the use of xlink: without them, we would have to repeat descriptions from a previous GML fragment
        inside the following <gml:operationParameter> elements.
      -->
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:value uom="urn:ogc:def:uom:EPSG::9102">28.583333</gml:value>
          <gml:operationParameter xlink:href="#epsg-parameter-8834"/> <!-- Latitude of topocentric origin -->
        </gml:ParameterValue>
      </gml:parameterValue>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:value uom="urn:ogc:def:uom:EPSG::9102">-80.583056</gml:value>
          <gml:operationParameter xlink:href="#epsg-parameter-8835"/> <!-- Longitude of topocentric origin -->
        </gml:ParameterValue>
      </gml:parameterValue>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:value uom="urn:ogc:def:uom:EPSG::9001">100</gml:value>
          <gml:operationParameter xlink:href="#epsg-parameter-8836"/> <!-- Ellipsoidal height of topocentric origin -->
        </gml:ParameterValue>
      </gml:parameterValue>
    </gml:Conversion>
  </gml:conversion>
  <!--
    The following block contains the Coordinate Reference System (CRS) associated to the observatory platform.
    In the GeoPose specification, this is called "outer frame." GeoPose restricts the CRS of outer frame
    to EPSG:4979, which is the CRS used in this example. But in GML any geodetic CRS could be used below.
    The XML file on the GitHub repository contains a copy of the full definition, but for this report we make
    the example shorter by using a link to EPSG registry.
  -->
  <gml:baseCRS xlink:href="https://epsg.org/api/v1/CoordRefSystem/4979/export?format=gml"/>
  <!--
    The following block defines the axes of the observatory platform as a Cartesian coordinate system.
    The platform origin is located at 28.583333°N, 80.583056°W and 100 meters above ellipsoid.
    Those values are specified in the <gml:conversion> block above.
  -->
  <gml:derivedCRSType codeSpace="EPSG">engineering</gml:derivedCRSType>
  <gml:coordinateSystem>
    <gml:CartesianCS gml:id="ObservatoryPlatformCS">
      <gml:identifier codeSpace="TB18-D025">ObservatoryPlatformCS</gml:identifier>
      <gml:name>Topocentric easting and northing</gml:name>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="PlatformAxisEast">
          <gml:identifier codeSpace="TB18-D025">PlatformAxisEast</gml:identifier>
          <gml:name>Topocentric East</gml:name>
          <gml:axisAbbrev>U</gml:axisAbbrev>
          <gml:axisDirection codeSpace="ISO">east</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:axis>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="PlatformAxisNorth">
          <gml:identifier codeSpace="TB18-D025">PlatformAxisNorth</gml:identifier>
          <gml:name>Topocentric North</gml:name>
          <gml:axisAbbrev>V</gml:axisAbbrev>
          <gml:axisDirection codeSpace="ISO">north</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:axis>
      <gml:axis>
        <gml:CoordinateSystemAxis uom="urn:ogc:def:uom:EPSG::9001" gml:id="PlatformAxisHeight">
          <gml:identifier codeSpace="TB18-D025">PlatformAxisHeight</gml:identifier>
          <gml:name>Topocentric height</gml:name>
          <gml:axisAbbrev>W</gml:axisAbbrev>
          <gml:axisDirection codeSpace="ISO">up</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:axis>
    </gml:CartesianCS>
  </gml:coordinateSystem>
</gml:DerivedCRS>

Figure 9 — Derived CRS for observatory platform (GML)

5.4.  Earth-centered inertial CRS

An Earth-centered inertial (ECI) coordinate reference system has its origin at the center of mass of the Earth and its axes oriented in directions fixed with respect to the stars. This contrasts with the “Earth-centered — Earth-fixed” (ECEF) frames used in classical geodesy. In geodesy, frames remain fixed with respect to Earth’s surface in its rotation, and then rotated with respect to stars. ECEF coordinates (typically latitude, longitude, and altitude) are convenient for representing the positions and velocities of terrestrial objects. But for objects in space, the equations of motion that describe orbital motion are simpler in a non-rotating frame such as ECI.

The extent to which an ECI frame is inertial is limited by the non-uniformity of the surrounding gravitational field. For example, the Moon’s gravitational influence on a high-Earth orbiting satellite is significantly different than its influence on Earth. Therefore observers in an ECI frame would have to account for this difference in acceleration in their laws of motion. The closer the observed object is to the ECI-origin, the less significant the effect of the gravitational disparity is.

Defining a truly inertial CRS is difficult. Only “quasi-inertial” approximations are used in practice. The 3D+ Standards Framework Engineering Report (OGC 22-036) discusses the different kinds of quasi-inertial CRS. For the purpose of this ER (Clause 5.4.1, Annex A.3.2, and VoyagerToObservatory.xml GML), this complexity is simplified with an hypothetical quasi-inertial CRS defined in such a way that the geodetic CRS is rotating at a constant speed relative to that quasi-inertial CRS.

The WKT below defines a quasi-inertial CRS (ECI) without specifying its relationship with a geodetic CRS (ECEF). The relationship between ECEF and ECI is not specified here because those two CRSs are defined by separated sets of observations (national survey versus astronomical observations). ECI is not defined in terms of ECEF, or conversely, ECEF is not defined in terms of ECI. Consequently, saying that one CRS is derived from the other is not appropriate. As such the use of DerivedCRS as in Clause 5.3 would not be correct. The relationship between ECI and ECEF is defined separately in Clause 5.4.1.

ENGINEERINGCRS["Quasi-inertial (Equatorial CS)",
  ENGINEERINGDATUM["ECI frame",
    ANCHOR["Earth center of mass"]],
  CS[Spherical, 3],
  AXIS["Declination (DEC)", north, ORDER[1], ANGLEUNIT["degree", 0.0174532925199433]],
  AXIS["Right Ascension (RA)", east, ORDER[2], ANGLEUNIT["degree", 0.0174532925199433]],
  AXIS["Radius (R)", up, ORDER[3], LENGTHUNIT["metre", 1]],
  REMARK["North stands for ""north celestial pole"" and East for ""6 hours from vernal equinox""."]]

Figure 10 — Quasi-inertial CRS (WKT)

COORDINATEOPERATION["Earth-Centered Inertial (ECI) to Earth Fixed (ECEF)",
  SOURCECRS[WKT from Clause 5.4 omitted for brevity],
  TARGETCRS[
    GEODETICCRS["Geocentric (Spherical CS)",
      DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84", 6378137, 298.257223563]],
      CS[Spherical, 3],
      AXIS["Latitude (P)", north, ORDER[1], ANGLEUNIT["degree", 0.0174532925199433]],
      AXIS["Longitude (L)", east, ORDER[2], ANGLEUNIT["degree", 0.0174532925199433]],
      AXIS["Radius (R)", up, ORDER[3], LENGTHUNIT["metre", 1]]]],
  METHOD["Time-dependent longitude rotation (equatorial CS)"],
  PARAMETER["Longitude axis rotation", 0, UNIT["degree", 0.0174532925199433]],
  PARAMETER["Rate of change of longitude axis rotation", 0.985609112046479, UNIT["degree/day", 1507.9644737231]],
  PARAMETER["Parameter reference epoch", 2000, UNIT["year", 31556925.445]]]

5.5.  Addition of temporal axis

All coordinate reference systems in the above clauses were spatial. ISO 19111 has a mechanism for adding a temporal axis to almost any CRS, but it assumes that the temporal dimension is independent from all other dimensions. This is not true when Einstein’s relativity becomes important, but that case will be discussed in Clause 8.5. This section assumes that relativistic effects can be ignored, in which case the process is as below:

• defines a TemporalCRS as a one-dimensional CRS, independent from the spatial CRS; and

• defines a CompoundCRS containing two components: the spatial CRS and above-cited temporal CRS.

The number of dimensions of the CompoundCRS is the sum of the number of dimensions of all its components. In this example the result is a four-dimensional CRS, but more components such as ParametricCRS could be added. The following WKT snippet shows a compound CRS using truncated Julian days for the temporal component. Truncated Julian is the time measured as days since May 24, 1968 at 00:00 UTC and is defined relative to Julian day (JD) as TJD = JD − 2440000.5. This epoch was introduced by NASA for the space program.

COMPOUNDCRS["Geocentric (Spherical CS) with time",
  GEODETICCRS[WKT from Clause 5.4.1 omitted for brevity],
  TIMECRS["Truncated Julian days",
    TDATUM["Truncated Julian epoch", TIMEORIGIN[1968-05-24T00:00:00Z]],
    CS[TemporalCount, 1],
    AXIS["Time", future, TIMEUNIT["day", 86400]]]]

Figure 12 — Compound CRS with a temporal component (WKT)

The addition of a temporal axis to an existing CRS is a complication when those CRSs need to be used with coordinate operations designed for spatial CRSs (for example map projections). This problem is handled by PassThroughOperation (Clause 5.6.3), which is defined in ISO 19111 and GML but not in WKT. Consequently, the WKT format is sufficient for expressing CRS with time axis but may not be sufficient for expressing coordinate operations with a temporal component.

COORDINATEMETADATA[
  ENGINEERINGCRS["A spacecraft-centred CRS for Voyager",
    See Clause 5.2 for the remaining of this WKT],
  EPOCH[2022.47]]

5.6.  Coordinate operations

ISO 19111:2019 distinguishes two kinds of coordinate operations (ignoring dynamic reference frame for now): conversions and transformations. A conversion can include translation, rotation, change of units, etc., but with a result which is still associated to the same reference frame, for example the same spacecraft. By contrast a transformation involves a change of reference frame, for example from one spacecraft to another one.

The distinction between conversion and transformation types allows type safety in DerivedCRS definitions. A conversion (but not a transformation) can be embedded in the definition of a derived CRS. That conversion is theoretically of infinite precision (ignoring rounding errors caused by floating-point arithmetic and numerical approximations caused by Taylor series expansions) because the derived CRS is by definition the result of applying that conversion on the base CRS. From the above, at most one conversion (ignoring differences in metadata) can exist between any given pair of CRSs in a self-consistent geodetic dataset. A derived CRS example was given in Clause 5.3.1.

By contrast, transformations between two CRSs are defined externally (for example in a database) and are not part of the CRS definition (ignoring BoundCRS, which is a compromise for a common but non-recommended practice). Transformations are determined empirically, for example by astronomical observations, and consequently have stochastic errors. Many transformations may exist between the same pair of CRSs, with different properties such as accuracy, domain of validity, etc.

In the context of Testbed 18, the operation for changing a coordinate tuple relative to a spacecraft into a coordinate tuple relative to another spacecraft (for example) would be a coordinate transformation, and the transformation parameters would be provided outside the CRS definitions, in an EPSG-like database. For map projection libraries using the EPSG database in a “late-binding” way, the existing software code should be applicable with few changes (Clause 9).

The coordinate operation type (conversion versus transformation) is not specified explicitly in the WKT format. The operation type is rather inferred from whether source and target datum are the same or different. The following example specifies a relationship between a pose and the extremity of a robotic arm. This example uses the quaternions parameters derived from the draft GeoPose standard as shown in Annex A.3.1. This example uses fixed parameter values, but different approaches for time-varying operations are introduced in Clause 5.4.1, Clause 5.6.1.2 and Clause 5.7.

COORDINATEOPERATION["Pose to extremity of robotic arm",
  SOURCECRS[GEOGRAPHICCRS[... complete definition omitted for brevity ...]],
  TARGETCRS[ENGINEERINGCRS[... complete definition omitted for brevity ...]],
  METHOD["Pose location and orientation by quaternion"],
  PARAMETER["Latitude of pose", 28.583333, ANGLEUNIT["degree", 0.0174532925199433]],
  PARAMETER["Longitude of pose", -80.583056, ANGLEUNIT["degree", 0.0174532925199433]],
  PARAMETER["Ellipsoidal height of pose", 100, LENGTHUNIT["metre", 1]],
  PARAMETER["Quaternion X",  0.20056154657066608, SCALEUNIT["unity", 1]],
  PARAMETER["Quaternion Y", -0.08111602541464237, SCALEUNIT["unity", 1]],
  PARAMETER["Quaternion Z",  0.36606032744426537, SCALEUNIT["unity", 1]],
  PARAMETER["Quaternion W", -0.90509396922613010, SCALEUNIT["unity", 1]],
  OPERATIONACCURACY[0.1]]

Figure 14 — Coordinate transformation from a pose to the extremity of a robotic arm

The above WKT does not define a CRS. It defines a relationship between two coordinate reference systems. An arbitrary number of relationships can exist and can be stored in an EPSG-like database. When a coordinate operation between the two CRSs is requested, map projection libraries such as PROJ or Apache SIS iterate over all coordinate operations found in the database for that pair of CRSs. Then an operation is selected using some heuristic rules, for example:

• preference given to the largest intersection between the temporal domain of validity of the coordinate operation and the time of interest specified by the user;

• preference given to the largest intersection between the geographic domain of validity of the coordinate operation and the area of interest specified by the user; and

• preference given to the operation having the smallest OPERATIONACCURACY value.

<gml:OperationMethod gml:id="PoseIrregularSeries">
  <gml:identifier codeSpace="TB18-D025">PoseIrregularSeries</gml:identifier>
  <gml:name>Pose irregular series</gml:name>
  <gml:formula>See TB18 D025.</gml:formula>
  <gml:sourceDimensions>3</gml:sourceDimensions>
  <gml:targetDimensions>3</gml:targetDimensions>
  <gml:parameter>
    <gml:OperationParameter gml:id="start-instant">
      <gml:identifier codeSpace="TB18-D025">:start-instant</gml:identifier>
      <gml:name>Start instant</gml:name>
    </gml:OperationParameter>
  </gml:parameter>
  <gml:parameter>
    <gml:OperationParameter gml:id="end-instant">
      <gml:identifier codeSpace="TB18-D025">:end-instant</gml:identifier>
      <gml:name>End instant</gml:name>
    </gml:OperationParameter>
  </gml:parameter>
  <!--
    Above "Start instant" and "End instant" parameters will have exactly one value.
    The parameters defined below are inside a group which can be repeated as many
    times as desired with different parameter values.
  -->
  <gml:parameter>
    <gml:OperationParameterGroup gml:id="pose-irregular-frame">
      <gml:identifier codeSpace="TB18-D025">pose-irregular-frame</gml:identifier>
      <gml:name>Frame specification (irregular)</gml:name>
      <gml:minimumOccurs>1</gml:minimumOccurs>
      <gml:maximumOccurs>1000</gml:maximumOccurs> <!-- Arbitrary limit in number of steps. -->
      <gml:parameter>
        <gml:OperationParameter gml:id="valid-time">
          <gml:identifier codeSpace="TB18-D025">valid-time</gml:identifier>
          <gml:name>Valid time</gml:name> <!-- As ISO-8601 or duration since Unix epoch. -->
        </gml:OperationParameter>
      </gml:parameter>
      <gml:parameter>
        <gml:OperationParameter gml:id="pose-location">
          <gml:identifier codeSpace="TB18-D025">pose-location</gml:identifier>
          <gml:name>Pose location</gml:name> <!-- As a sequence of 3 coordinate values. -->
        </gml:OperationParameter>
      </gml:parameter>
      <gml:parameter>
        <gml:OperationParameter gml:id="yaw">
          <gml:identifier codeSpace="TB18-D025">yaw</gml:identifier>
          <gml:name>Yaw</gml:name>
        </gml:OperationParameter>
      </gml:parameter>
      <gml:parameter>
        <gml:OperationParameter gml:id="pitch">
          <gml:identifier codeSpace="TB18-D025">pitch</gml:identifier>
          <gml:name>Pitch</gml:name>
        </gml:OperationParameter>
      </gml:parameter>
      <gml:parameter>
        <gml:OperationParameter gml:id="roll">
          <gml:identifier codeSpace="TB18-D025">roll</gml:identifier>
          <gml:name>Roll</gml:name>
        </gml:OperationParameter>
      </gml:parameter>
    </gml:OperationParameterGroup>
  </gml:parameter>
</gml:OperationMethod>

<gml:Conversion gml:id = "PoseToRoboticArm">
  <gml:identifier codeSpace="TB18-D025">PoseToRoboticArm</gml:identifier>
  <gml:name>Pose to extremity of a robotic arm</gml:name>
  <gml:scope>TestBed 18 demonstration.</gml:scope>
  <gml:method xlink:href="#PoseIrregularSeries"/> <!-- Reference to above <gml:OperationMethod>. -->
  <gml:parameterValue>
    <gml:ParameterValue>
      <gml:stringValue>2022-10-20T06:00:00Z</gml:stringValue>
      <gml:operationParameter xlink:href="#start-instant"/>
    </gml:ParameterValue>
  </gml:parameterValue>
  <gml:parameterValue>
    <gml:ParameterValue>
      <gml:stringValue>2022-10-21T18:00:00Z</gml:stringValue>
      <gml:operationParameter xlink:href="#end-instant"/>
    </gml:ParameterValue>
  </gml:parameterValue>
  <gml:parameterValue>
    <!--
      Yaw/pitch/roll parameters for the first time step.
    -->
    <gml:ParameterValueGroup>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:stringValue>2022-10-20T06:00:00Z</gml:stringValue>
          <gml:operationParameter xlink:href="#valid-time"/>
        </gml:ParameterValue>
      </gml:parameterValue>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:valueList uom="urn:ogc:def:uom:EPSG::9001">4 6 2</gml:valueList> <!-- UoM is metres. Values are dummy. -->
          <gml:operationParameter xlink:href="#pose-location"/>
        </gml:ParameterValue>
      </gml:parameterValue>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:value uom="urn:ogc:def:uom:EPSG::9122">40</gml:value> <!-- UoM is degrees. Value is dummy. -->
          <gml:operationParameter xlink:href="#yaw"/>
        </gml:ParameterValue>
      </gml:parameterValue>
      <!-- Pitch and Roll parameters omitted for brevity. -->
      <gml:group xlink:href="#pose-irregular-frame"/>
    </gml:ParameterValueGroup>
  </gml:parameterValue>
  <gml:parameterValue>
    <!--
      Yaw/pitch/roll parameters for the second time step.
      This is the same group than above, repeated with different parameter values.
    -->
    <gml:ParameterValueGroup>
      <gml:parameterValue>
        <gml:ParameterValue>
          <gml:stringValue>2022-10-20T17:00:00Z</gml:stringValue>
          <gml:operationParameter xlink:href="#valid-time"/>
        </gml:ParameterValue>
      </gml:parameterValue>
      <!-- Pose location, Yaw, Pitch and Roll parameters omitted for brevity. -->
      <gml:group xlink:href="#pose-irregular-frame"/>
    </gml:ParameterValueGroup>
  </gml:parameterValue>
  <!--
    More time steps could be added by continuing to repeat <gml:ParameterValueGroup>.
  -->
</gml:Conversion>

CONCATENATEDOPERATION["Voyager spacecraft to Observatory platform",
  SOURCECRS[ENGINEERINGCRS["A spacecraft-centred CRS for Voyager", …]],
  TARGETCRS[GEODETICCRS["Local Tangent Plane - East North Up (LTP-ENU)", …]],
  STEP[COORDINATEOPERATION[
    SOURCECRS[ENGINEERINGCRS["A spacecraft-centred CRS for Voyager", …]],
    TARGETCRS[INERTIALCRS["Ecliptic heliocentric CRS (Spherical)", …]],
    METHOD["Trajectory to Conventional frame"], …]],
  STEP[COORDINATEOPERATION[
    SOURCECRS[INERTIALCRS["Ecliptic heliocentric CRS (Spherical)", …]],
    TARGETCRS[INERTIALCRS["Quasi-inertial (Equatorial CS)", …]],
    METHOD["To circular orbit (Spherical domain)"], …]],
  STEP[COORDINATEOPERATION[
    SOURCECRS[INERTIALCRS["Quasi-inertial (Equatorial CS)", …]],
    TARGETCRS[GEODETICCRS["WGS 84 (Spherical)", …]],
    METHOD["Longitude axis rotation"], …]],
  STEP[COORDINATEOPERATION[
    SOURCECRS[GEODETICCRS["WGS 84 (Spherical)", …]],
    TARGETCRS[GEOGRAPHICCRS["WGS 84", …]],
    METHOD["Geographic/spherical conversions"], …]],
  STEP[COORDINATEOPERATION[
    SOURCECRS[GEOGRAPHICCRS["WGS 84", …]],
    TARGETCRS[GEODETICCRS["Local Tangent Plane - East North Up (LTP-ENU)", …]],
    METHOD["Geographic/topocentric conversions"], …]]

<gml:PassThroughOperation gml:id="SphericalToGeographic4D">
  <gml:identifier codeSpace="TB18-D025">SphericalToGeographic4D</gml:identifier>
  <gml:name>Earth Fixed (ECEF) spherical to geographic with time</gml:name>
  <gml:scope>TestBed 18 demonstration.</gml:scope>
  <gml:sourceCRS xlink:href="#GeodeticAndTimeCRS"/>
  <gml:targetCRS xlink:href="#GeographicAndTimeCRS"/>
  <gml:modifiedCoordinate>1</gml:modifiedCoordinate> <!-- Index values start at 1. -->
  <gml:modifiedCoordinate>2</gml:modifiedCoordinate>
  <gml:modifiedCoordinate>3</gml:modifiedCoordinate>
  <gml:coordOperation>
    <gml:Conversion gml:id="SphericalToGeographic">
      <gml:identifier codeSpace="TB18-D025">SphericalToGeographic</gml:identifier>
      <gml:name>Earth Fixed (ECEF) spherical to geographic</gml:name>
      <gml:scope>TestBed 18 demonstration.</gml:scope>
      <gml:method>
        <gml:OperationMethod gml:id="GeographicSphericalConversion">
          <gml:identifier codeSpace="TB18-D025">GeographicSphericalConversion</gml:identifier>
          <gml:name>Geographic/spherical conversions</gml:name>
          <gml:formula>See TB18 D025.</gml:formula>
          <gml:sourceDimensions>3</gml:sourceDimensions>
          <gml:targetDimensions>3</gml:targetDimensions>
        </gml:OperationMethod>
      </gml:method>
    </gml:Conversion>
  </gml:coordOperation>
</gml:PassThroughOperation>

5.7.  Spacecraft trajectory

The position of a spacecraft may be too complex for a description using simple formulas with a few parameters. Instead, providing the spacecraft trajectory in a file may be more convenient. This approach is similar to the use of datum shift grids (NADCON, NTv2, etc.) but with vector data instead of gridded data. In Testbed 18, the participants invented a new operation method named “Trajectory to Conventional frame” with a single parameter, the name of a file encoded in Moving Feature CSV format (OGC 14-084r2). The following constraints shall hold:

• operation source CRS = CRS attached to the moving feature (Voyager in our scenario); and

• operation target CRS = CRS of the Moving Feature CSV file (Heliocentric in our scenario).

Spacecraft yaw, roll, and pitch can be provided as attributes of the moving features. Using the spacecraft location relative to the Sun together with yaw/roll/pitch information, the coordinates of an object relative to Voyager can be transformed to heliocentric coordinates.

Below is an example of moving features CSV files. The first line declares the CRS, bounding boxes (in CRS units), start time, and end time. The second line declares the columns, with yaw/roll/pitch attributes as the last three columns (the corresponding values in this example are random). The HeliocentricCRS identifier can be used for looking up the complete CRS definition in the VoyagerToObservatory.xml file. CRS axes are Ecliptic latitude (degrees), Ecliptic longitude (degrees), and Distance from Sun (km). This example contains only one feature — the Voyager spacecraft — but an arbitrary number of features could be stored in the same file.

@stboundedby, urn:ogc:def:crs:TB18-D025::HeliocentricCRS, 3D, -37.524474 289.911611 19695916300, -37.524474 289.911611 19695917500, 2022-10-31T18:52:00Z, 2022-10-31T19:00:00Z, absolute
@columns, mfidref, trajectory, yaw,xsd:decimal, pitch,xsd:decimal, roll,xsd:decimal
Voyager, 2022-10-31T18:53:00Z, 2022-10-31T18:54:00Z, -37.524474 289.911611 19695916300, 10, 5, 20
Voyager, 2022-10-31T18:54:00Z, 2022-10-31T18:55:00Z, -37.524474 289.911611 19695916800, 10, 4, 21
Voyager, 2022-10-31T18:55:00Z, 2022-10-31T18:56:00Z, -37.524474 289.911611 19695917500, 10, 5, 20

Figure 19 — Spacecraft coordinates (Moving Features CSV)

The demonstration Java application (Clause 9.1) defines a “Trajectory to Conventional frame” operation method capable to read and use this CSV file. The coordinate operation is implemented by the TrajectoryToConventionalFrame.java file and added to the map projection library using the plugin mechanism described in Clause 8.9.

5.8.  Relativist CRS

The proposals in all previous sections (except Clause 5.7.2) ignored Einstein’s relativity. In all previous CRS definitions, the temporal dimension was independent from the spatial dimensions. Those dimensions were defined in separated SingleCRS components, then those components were assembled in a CompoundCRS. In the ISO 19111 model, CompoundCRS is not associated with a coordinate system (Figure 1). Only the CRS components (the spatial and temporal CRSs) are associated with a coordinate system. Consequently, there is no 4-dimensional coordinate system. Only a 3-dimensional CartesianCS (or other type) followed by a 1-dimensional TemporalCRS.

For objects moving at sufficient speeds relative to an observer, the above specified separation is not desirable. The rate of change in the three spatial dimensions (the velocity) has an impact on the rate of change in the fourth dimension (the frequency of clock ticks). For doing coordinate transformations in the context of Einstein’s relativity, the two coordinate systems (3D + 1D) need to be replaced by a single coordinate system (4D). A convenient 4D coordinate system for relativist calculations is the Minkowski space. However, there is no “Minkowski coordinate system” in the current ISO 19111. Adding that coordinate system would be an extension, which is proposed in Clause 8.5.

Given two coordinate reference systems named “A” and “B,” if A and B are in different gravity fields and/or are moving fast relative to each other, then (x, y, z, t) coordinate tuples can be transformed from A to B in two different ways:

By formula: If the “gravity terrain” is simple enough, an application could do the mathematics by itself. A and B would need to be SingleCRS instances associated to a Minkowski coordinate system. A new OperationMethod instance would need to be defined in a way similar to Clause 5.3.1, then used in the definition of a coordinate operation similar to Clause 5.6.

By gridded data: If the “gravity terrain” is too complex, then some data producer would need to run a model and save the output in a “datum shift grid.” This is similar to NADCON grids or NTv2, but in 4 dimensions instead of 2. The Gridded Geodetic data eXchange Format (GGXF) could be used (Clause 5.7.2).

5.9.  Spacecraft deformations

The OGC Temporal Domain Working Group (DWG) has begun development of a discussion paper on dynamic features. Dynamic features are moving features that can change location, shape, and state as a function of time. A complete description of the spacecraft shape and state is out of scope for this Engineering Report. However, the consequence of those changes on coordinate values is considered. For example, if an EngineeringCRS is attached to the spacecraft (Clause 5.2) and if all measurements on the spacecraft surface are relative to an anchor point, then the coordinates of another fixed point may change with time because, for example, of thermal dilatation. In this scenario, the change in coordinate values is not caused by an “active” displacement of a point. The location is still the same point on the spacecraft surface. The coordinate change is rather caused by a deformation of the support. In ISO 19111, this is called “point motion operation” and is a third kind of operation in addition to conversion and transformation (Clause 5.6). A reference frame subject to point motions is called a dynamic reference frame. By contrast, a reference frame which is not dynamic is called a static reference frame.

The difference between a dynamic reference frame and a dynamic feature is that the deformations of a dynamic feature can be expressed in a static reference frame, which is external to the feature. For example, the changing shape of a hurricane may be described in a static geographic CRS. A reference frame can become dynamic when it is attached to the shape of the dynamic feature.

Point motion operations are supported by the GGXF format (Clause 5.7.1). While it was designed primarily for tectonic plate movements, handling thermal dilatations of spacecraft should be possible as a problem similar to tectonic plate movements on Earth.

6.1.  Outer frame

GeoPose uses an outer frame as the starting point of a frame graph. In GeoPose 1.0, that outer frame is fixed to a three-dimensional “World Geodetic System 1984” CRS with latitude, longitude, and ellipsoidal height axes (EPSG:4979). In the ISO 19111 model, there is no distinction between outer and inner frames. Any frame can be associated to the source CRS of a coordinate operation. Consequently, the restriction to EPSG:4979 does not apply.

GeoPose allows the definition of a Local Tangent Plane (LTP), which is tangential to the geoid or ellipsoid at a given latitude and longitude. This plane can be oriented in different ways such as (yaw, pitch, roll) angles or by quaternion, and the orientation can be time dependent. If defined, this tangent plane is an inner frame derived from the outer frame in the frame graph. However, as previously stated, in ISO 19111 there is no concept of inner and outer frames. Instead, the important criterion is whether the “inner frame” is referenced to the same object than the “outer frame”. For example, if the “inner frame” is a land area (or maybe a non-moving platform fixed to Earth), then the coordinate system is changed for convenience but the CRS is still associated (even if indirectly) to the same datum, for example WGS84. In such case there is a derived CRS as defined in Clause 5.3. Otherwise there is a coordinate transformation as defined later in this section.

If a local tangent plane is defined by a DerivedCRS, then there does not need to be an explicit part of the frame graph. Since ISO 19111 has no restriction about the source CRS of a coordinate operation, a DerivedCRS can be the “outer frame” of a frame graph. Derived CRSs and their relationship with their base CRSs have well-defined encoding in WKT (Clause 5.3) and GML.

Another way to decide whether to use a DerivedCRS is to look how its relationship with the base CRS was defined. As a rule of thumb (not an absolute criterion), if the relationship has been determined by observations, is subject to stochastic errors, and can be revised at any time, then there probably is a coordinate transformation. Conversely if the relationship has been established by definition (potentially after observation, before an authority decided to fix some parameter values as the official ones), then there probably is a coordinate conversion. In the former case (transformation), DerivedCRS cannot be used. In the latter case (conversion), a DerivedCRS can be used, but it still optional. In any case, the distinction between conversion and transformation carries semantic information which is absent from GeoPose.

6.2.  Inner frame

An inner frame is defined by the coordinate operation from another (possibly outer) frame to the inner frame. This is similar to a ISO 19111 derived CRS, except that the coordinate operation can be a transformation. GeoPose defines various classes (YawPitchRollAngles, AtAndUpVector, etc.) for storing the parameters required for the definition of coordinate operations. Each class is specialized for some specific parameters, for example rotation angles. This is in contrast with ISO 19111, which has an approach more similar to ISO 19109 FeatureType. In ISO 19109, a FeatureType is a meta-class: a class used for defining other classes. This is similar to reflection in the Java programming language, or to the elements used for defining a schema in XML or JSON. In ISO 19111, the OperationParameter and OperationParameterGroup types (Figure 7) act like meta-classes: they define the set of expected parameters for a given coordinate operation, much as a JSON schema defines the set of properties for a given JSON element. When applied to strongly typed programming languages such as Java, meta-classes are flexible because they do not require the definition of new classes in the programming language for each new kind of operation. Instead, new operations are described by instances of ISO 19111 OperationMethod and OperationParameter types, which in turn can be materialized by a few rows in the EPSG database. Those definitions need to be provided only once and can be referenced by xlink in a GML document or by authority code in a WKT definition.

GeoPose defines many sets of parameters, from basic to more advanced use cases. But they are all different ways to specify translations and rotations. This is similar to EPSG providing different ways for defining the same map projection. For example, EPSG provides three different ways to define a Mercator projection:

• Mercator (variant A) which accepts a Scale factor at natural origin parameter (among others);

• Mercator (variant B) which accepts a Latitude of 1st standard parallel parameter (among others); and

• Mercator (variant C) which accepts Latitude of 1st standard parallel and Latitude of false origin parameters (among others).

All those variants are only different ways to process the parameters. After parameter processing, those variants can be implemented by the same code. The same remark applies to operations defined by GeoPose. This is demonstrated in Annex A.3.1, which lists GeoPose operations reworded as definitions using ISO 19111 objects. The following GeoPose operation methods are defined:

• Basic-YPR (yaw, pitch, roll) ⟶ “Pose rotation by YPR;” and

• Basic-Quaternion ⟶ “Pose location and orientation by quaternion.”

The Advanced GeoPose method is excluded because it can be done by using a DerivedCRS as the outer frame (Clause 6.1), and then using the Basic-Quaternion operation method (Clause 5.6). The GEODETICCRS["Local Tangent Plane - East North Up (LTP-ENU)"] element in Clause 5.3 corresponds to a first inner frame in GeoPose. Then the COORDINATEOPERATION["Pose to extremity of robotic arm"] element in Clause 5.6 corresponds to a second inner frame in GeoPose.

6.3.  Chain and frame graph

A chain is an outer frame and a sequence of transformations to a final innermost frame. A frame graph is a structured modeling of the pose relationship between frames (nodes) and transforms (edges) in a graph structure. Those two constructs can be modeled in ISO 19111 by CoordinateOperation and ConcatenatedOperation. See Clause 5.6.2 for an example.

6.4.  Time-varying parameters

GeoPose provides two ways to define a coordinate operation with time-dependent parameters. One way is termed “regular” and is built with a sequence of frames with a constant inter-pose duration. The second way is termed “irregular” and supports variable inter-pose durations. In both cases, there is a block of parameters which is repeated with different values for each GeoPose frame. In the ISO 19111 model, a repeating set of parameters can be represented by a OperationParameterGroup (Figure 7). A group can contain for example the “valid time”, “translation,” and “rotation” parameters. Parameter groups derived from the GeoPose standard are given in Annex A.2. Then a coordinate operation can be created with a set of parameters containing a OperationParameterGroup instance for each valid time. Clause 5.6.1.2 shows an example in GML (the parameter group construct is not available in WKT) using the “Pose irregular series” transformation. A complete GML document is provided by the PoseToRoboticArm.xml file in the Testbed 18 GitHub repository.

The following are two differences in the ISO 19111 model compared to the GeoPose model.

• Translation and rotation are not parts of the frame definition. They are parts of the coordinate operation between a pair of frames.

• The set of parameters for all valid times are grouped in a single operation. The above example does not define a new frame for each valid time.

Note that all examples shown in this document up to this point use existing ISO 19111 and GML standards. No extension to existing standards has been required yet.

6.5.  GeoPose encoding of CRS

GeoPose encodes coordinate reference systems (CRS) in JSON documents with a small block containing an authority, an identifier, and a list of parameters. This block and its properties are specific to the GeoPose draft standard. An issue with the current GeoPose approach is that it is implementation dependent. The authority property can be interpreted as the language in which parameters are encoded. If the authority is PROJ, then the parameters are a PROJ string and only PROJ (or other software having some compatibility layer) can understand those parameters. The same issue exists with EPSG codes as well. If a software does not have a connection or a copy of EPSG geodetic dataset, it cannot understand the CRS.

The GeoPose structure could be replaced by a ISO 19111 structure encoded in JSON. There is no OGC standard yet for JSON encoding of ISO 19111, but a proposal based on PROJ-JSON is on the table. The ISO 19111 structure contains all information currently found in GeoPose structure, plus an implementation neutral description of the CRS. For example, no matter the authority, all CRSs have a set of axes. The ISO 19111 structures enumerate those axes together with other information such as datum. Those structures do not necessarily contain all metadata provided by the authority, but they contain sufficient information for allowing coordinate operations even if the software does not know the authority. CRS encoding derived from ISO 19111 such as GML, WKT, or PROJ-JSON make possible to:

• use the authority code if the software understands the authority;

• or otherwise fallback on information embedded in the CRS definition.

7.1.  Getting matrix values

If the coordinate operation is an affine transform, then the Jacobian matrix is the same everywhere. But for any non-linear operation, the Jacobian matrix can be different at every spatiotemporal location. Consequently, Jacobian matrices must be computed in a way similar to the computations done for transforming positions. ISO 19111 currently defines the following method on CoordinateOperation:

CoordinateSet transform(CoordinateSet positions);

In addition to above method, the following method would also be needed. This method already exists in GeoAPI MathTransform interface, inherited from OGC 01-009:

Matrix derivative(DirectPosition position)

Those two methods exist for ISO 19111 and OGC 01-009, respectively. A possible extension to those APIs is presented in Clause 8.6.2.

7.2.  Example: transforming a BBOX

Some numerical problems, such as the search for minimum value, are solved by iterative computations of a function. Often, the search converges faster when the function derivative is known (e.g., Gradient descent algorithm). This is also true for the transformation of a bounding box from one CRS to another CRS. The algorithm needs to find the minimum and maximum values of the transformed shape. If the transformation is non-linear, then those minimums and maximums are not necessarily the corner coordinates. For example, the blue shape in the figure below is the result of transforming a bounding box. As one can see, the lowest part of the blue shape is lower than all corners. The transform algorithm needs to find the minimal y coordinate, which appears somewhere between the corners. A brute force approach would be to sample a large number of points between the corners and retain the smallest value, with the hope that the sampling is dense enough for including a point close to the real minimum. But a more efficient approach is to use the derivatives at the two lower corners for approximating a cubic polynomial (y = A + B⋅λ + C⋅λ² + D⋅λ³), illustrated by a red curve in the figure below.

Figure 21 — Use of derivatives for BBOX transformations

The cubic polynomial can easily be solved for finding the λ coordinate of the minimal y value according that approximation. That minimum computed form cubic approximation is usually located close to the real minimum of the transformed shape. If desired the Jacobian matrix can be evaluated recursively at the new λ coordinates for getting even closer. Then the real y value can be computed at the λ coordinate of the minimum.

The performance advantage of using Jacobian matrices for BBOX transformations increases quickly with the number of dimensions. A brute force approach sampling 100 points on each side would need to sample 4×100 points for a two-dimensional rectangle, 100×100×6 points for a three-dimensional cube, and millions of points for a four-dimensional hypercube (3D + time). In practice, the number of points can be greatly reduced when considering some dimensions as independent (e.g., time), or if the transformation is linear in at least some dimensions. But in a relativist system the 4 dimensions are not independent and the transformations are not always linear.

7.3.  Example: datum shift

Some transformations are too complex for expression by a mathematical formula. Instead, a grid is provided with displacement vectors from coordinates in the source CRS to coordinates in the target CRS. This grid can be encoded in a GGXF file (Clause 5.7.1), for example.

Figure 22 — Displacement vectors (source: GGXF group)

Transforming coordinates from a source CRS to a target CRS is a straightforward process as follows.

1. Get the displacement vectors at source coordinates.

2. Compute target coordinates = source coordinates + displacement vector.

However, the reverse operation (i.e., from target to source) is much more complicated. A displacement vector cannot be directly used for given target coordinates because the domain of the grid is in the source CRS, not the target CRS. The coordinates in source CRS (i.e., the location where to fetch a displacement vector) are not known since they are precisely the values we want to compute. For solving this problem, NADCON (North American Datum Conversion) uses an iterative approach based on the assumption that source coordinates are approximately equal to target coordinates, with only small displacement vectors between them. This assumption allows the following algorithm (using vector arithmetic).

1. Use target coordinates as initial interim source coordinates.

2. Get the displacement vector at interim source coordinates.

3. Compute interim source coordinates = target coordinates − displacement vector.

4. Repeat from step 2 until the change between two iterations is smaller than the desired accuracy.

However, the above algorithm does not converge if there are relatively large variations in vector magnitudes and directions. The interim source coordinates may move chaotically around the true source coordinates and can even go further and further at each iteration step. A more stable (but still not perfect) algorithm is as shown below.

1. Use target coordinates as initial interim source coordinates.

2. Get the displacement vector at interim source coordinates.

3. Compute interim source coordinates = target coordinates − displacement vector.

4. Get new displacement vector at new interim source coordinates.

5. Compute error in target = (new interim source coordinates + new displacement vector) − target coordinates.

6. Compute error in source = J⁻¹ × error in target where J is defined below.

7. Translate the interim source coordinates by a vector in the opposite direction of error in source.

8. Repeat from step 3 until the change between two iterations is smaller than the desired accuracy.

In a nutshell, compute at each iteration the error between the estimated position and the desired position. Then correct the position by moving in the direction opposite to the error. The problem is that the error is measured in target coordinates while we need to correct the source coordinates. Consequently, the error in target vector needs to be converted to an error in source vector. If the Jacobian matrix at the position of interim source coordinates is known, then the conversion can be done as below where (φ,λ) are source coordinates, (x, y) are target coordinates and ∆ are differences between estimated and actual values:

Note that if target coordinates are approximately equal to source coordinates, then the Jacobian matrix is close to an identity matrix. Consequently (∆φ, ∆λ) ≈ (∆x, ∆y) and step 7 in above algorithm simplifies as shown below.

• New interim source coordinates defined in step 7: (φ,λ) — (∆φ, ∆λ) ≈ (φ,λ) — (∆x, ∆y).

• From step 5: (∆x, ∆y) = (φ,λ) + (Dx, Dy) — (x, y) where D is the displacement vector.

• After substitution and simplification, interim source coordinates = (x, y) — (Dx, Dy).

The result of the above simplification is exactly step 3 in the NADCON algorithm. From the point of view of the algorithm using Jacobian matrix, datum shift algorithms such as NADCON have an implicit assumption that the Jacobian matrix is close to identity everywhere. This is true for most geodetic datum shifts, but it is not true for transformations based on grids having stronger curvatures or deformations.

In general relativity, gravity fields change the spatial and temporal coordinates all together. The gravity field may be too complex to be described by a few equations, as seen for example with Earth geoids. Given two reference frames in different gravity fields, transformations from one frame to the other frame could be described empirically by a grid similar to datum shift grids used in geodesy (Clause 5.7.2). If curvature is strong, classical algorithms used in geodesy may not converge. Jacobian matrices can resolve the convergence problem in some cases (not always, it also has its limits). Even when convergences are possible without them, as seen in previous section, Jacobian matrices can speed calculations by reducing the number of iterations, especially in spaces of more than two dimensions.

However, computing the Jacobian matrix has a cost. In a datum shift grid, it requires computing finite differences with neighbor nodes and the results are not exact representations of the derivatives at the desired location. Better accuracy, and possibly better performance, can be achieved if the Jacobian matrices are computed by the data producer and stored in the grid. The Gridded Geodetic data eXchange Format (GGXF), in development at OGC, is technically capable to store matrices.

7.4.  Lorentz transformations

In the context of special relativity, Lorentz transformations can be represented by a Jacobian matrix as discussed in Clause 5.8.1. In the context of general relativity, the Jacobian matrix is valid only locally and may change at every spatiotemporal position. But this is similar to the map projection case discussed in Clause 7, and the same API as Clause 7.1 (i.e., the derivative(…) method) can apply if the following conditions hold.

• The relative velocity between source CRS and target CRS is known (it may be a parameter of the coordinate operation).

• Coordinates and property values to transform are in units of the source CRS, not proper lengths (Clause 5.8.2), so the velocity of individual feature instances can be ignored.

The Jacobian matrix obtained under those conditions allows transformation of lengths observed in source CRSs to lengths observed in target CRSs. None of those lengths need to be proper lengths. If transformation from/to proper lengths is desired, then the transformation methods need an additional API discussed in Clause 8.6.2.

8.1.  Celestial body information

Technically, geodetic frames could be associated to any celestial body (if we ignore the etymology of “Geo” prefix in class names — see Clause 8.3). In practice, there is often an implicit assumption that geodetic frames are associated to the Earth. There is currently no way in the ISO 19111 standard to unambiguously identify the celestial body of a reference frame. A process could be to parse the reference frame name or the anchor, but such an approach is fragile because those attributes are free text strings. A more systematic approach would be to introduce a new class, for example named CelestialBody and add an association from GeodeticReferenceFrame to CelestialBody. Expressed in WKT format, this proposal would be the element in bold characters below:

GEODCRS["Mars (2015) / planetographic",
  DATUM["Mars (2015)",
    CELESTIALBODY["Mars"],
    ELLIPSOID["Mars", 3396190, 169.8944472236, LENGTHUNIT["metre", 1]],
    ANCHOR["Viking 1 lander : 47.95137 W"],
    PRIMEM["Airy", 0, ANGLEUNIT["degree", 0.017453292519943295]]],
  CS[ellipsoidal, 2],
    AXIS["Latitude (B)", north, ORDER[1]],
    AXIS["Longitude (L)", east, ORDER[2]],
    ANGLEUNIT["degree", 0.017453292519943295],
  SCOPE["For horizontal positioning on Mars."]]

Figure 23 — Celestial body proposal in a WKT definition

The PROJ 6+ library (Clause 9) already implements something similar. That library stores EPSG data in its own database schema, which is slightly different than the EPSG schema. The PROJ database schema contains a field for celestial body information.

In the ISO 19111 abstract model, the CelestialBody could be one of the following types.

• A subclass of ISO 19111 IdentifiedObject class. The advantage is to introduce no new dependency to another standard, which is inconvenient because the celestial body would be identified only by a single name, which may be insufficient.

• A subclass of ISO 19112 Location class (from Referencing by identifiers). The advantage is to enable the use of structured identifiers like addresses. For example, the celestial body could be identified by a “Galaxy → Star → Planet” hierarchy of names, however an inconvenience would be to introduce a dependency from ISO 19111 to ISO 19112.

Like ISO 19111, the terminology used in ISO 19112 assumes a location on Earth. The use of ISO 19112 depends on whether the use of GeoXXX objects is acceptable for bodies other than Earth. This issue is discussed in Clause 8.3.

8.2.  Non-geographic domain of validity

The BBOX element in the WKT format accepts only Earth latitudes and longitudes. Other ways to specify a domain of validity are needed, for example as a distance from an arbitrary celestial body. This is important for allowing map projection libraries to select the most appropriate coordinate operation for a pair of CRSs when more than one operation is available in the database (see Clause 5.6).

For simplicity reasons the latitude and longitude restriction is specific to the WKT format. The ISO 19111 abstract model does not use specifically the GeographicBoundingBox or GeographicExtent class. Instead, the model uses the more generic Extent class, which allows custom Extent subclasses. The steps for extending the domain of validity to spatial environments could be as follows.

• Create a new class as a specialization of ISO 19115 Extent.

• Create a new WKT keyword for representing the new extent class.

8.3.  Use of geocentric term

ISO 19111 uses “geographic” and “geocentric” terms in places where it could apply to any celestial bodies. For example, the AxisDirection code list contains geocentricX, geocentricY, and geocentricZ values. There are at least three possible ways to extend ISO 19111 axis directions to planetary CRS.

• Define new code list values such as heliocentricX, heliocentricY, and heliocentricZ for other astronomical bodies.

• Accept the use of geocentricX, geocentricY, and geocentricZ for any astronomical bodies, despite the “geo” semantic.

• Do a compromise with a single set of new axis directions for a whole category of astronomical bodies other than Earth. For example planetocentricX, planetocentricY and planetocentricZ.

The first alternative would make implementations more complicated because it would hide the fact that those axes are used with the same formulas (such as spherical to Cartesian coordinate system conversions) regardless of the astronomical body. The second alternative is the opposite extreme. The third alternative is between the two “extremes”.

The above discussion was only about one code list. But the “geo” prefix is spread across OGC/ISO standards. “Geo” appears in class names of ISO 19111 (GeographicCRS, GeodeticCRS, GeodeticReferenceFrame), of ISO 19115 (GeographicExtent, GeographicBoundingBox, GeographicDescription, GeolocationInformation, Georectified, Georeferenceable), and of ISO 19107 (Geodesic) to cite just a few. Earth-related prefixes also appear in attributes or associations (geographicIdentifier, alternativeGeographicIdentifier, geographicExtent, geographicElement, geolocationInformation, geoidModel, greenwichLongitude, etc). Duplicating all above-cited classes with different names, properties, and code list values would be significant work not only for the planetary working group, but also for the working groups of all impacted standards and for the implementers.

8.4.  Inertial CRS type

Clause 5.4 defined an Earth-centered quasi-inertial CRS using the EngineeringCRS type. This is because no other ISO 19111 type would have been acceptable. However, EngineeringCRS is still not quite right because this CRS is intended for local uses or vessels. A new CRS type is desirable, with a structure (class, associations, and attributes) similar to GeodeticCRS. The new CRS type could be named InertialCRS and its datum class could be named InertialDatum. Whether InertialDatum should contain an ellipsoid (the existing Ellipsoid class could be reused) is an open question. An ellipsoid may not be needed if the CRS accepts only three-dimensional Cartesian and spherical coordinate systems and is not intended for referencing objects on planet surface.

Using this InertialCRS proposal, the definition shown in Quasi-inertial CRS (WKT) figure of Clause 5.4 would become as below. This proposal also uses the CelestialBody proposal presented in Clause 8.1, and keep an ellipsoid inside the datum for this version.

INERTIALCRS["Quasi-inertial (Equatorial CS)",
  INERTIALDATUM["ECI frame",
    CELESTIALBODY["Earth"],
    ELLIPSOID["WGS 84", 6378137, 298.257223563]],
  CS[Spherical, 3],
  AXIS["Declination (DEC)", north, ORDER[1], ANGLEUNIT["degree", 0.0174532925199433]],
  AXIS["Right Ascension (RA)", east, ORDER[2], ANGLEUNIT["degree", 0.0174532925199433]],
  AXIS["Radius (R)", up, ORDER[3], LENGTHUNIT["metre", 1]],
  REMARK["North stands for ""north celestial pole"" and East for ""6 hours from vernal equinox""."]]

Figure 24 — Quasi-inertial CRS (proposed WKT extension)

8.5.  Minkowski coordinate system

Each Coordinate Reference System (CRS) is associated with a Coordinate System (CS). The subtypes defined by ISO 19111 such as CartesianCS and SphericalCS can be at most three-dimensional. The TemporalCS component, if desired, must be added separately as a dimension independent from the spatial ones (Clause 5.5).

For coordinate operations involving special or general relativity, the mathematical formulas require that space and time are handled together. To represent this requirement in the ISO 19111 model, a new CoordinateSystem subtype is required. The proposal is to define MinkowskiCS class as a coordinate system with the following constraints.

• The number of dimensions must be exactly 4.

• All dimensions must have the same unit of measurement (UoM). Temporal coordinates shall be converted to length unit by multiplying by the speed of light.

For each type of coordinate system, there is a different set of formulas for computing geometric properties such as distance. For example, the distance between two points is computed using the Pythagorean formula if the coordinate system is Cartesian, or a different formula if the coordinate system is polar. Like any coordinate system, the use of a Minkowski coordinate system also implies a particular formula for computing interval between two spacetime events:

(Δs)² = (Δct)² − (Δx)² − (Δy)² − (Δz)²

For the purpose of this Engineering Report, the key information to retain is that the use of MinkowskiCS implies not only a 4-dimensional space, but also homogeneity in axis units and a set of mathematical formulas that apply in this space.

8.6.  Relativist coordinate operations

Lorentz contractions (Clause 5.8.1) depend on relative velocity. Consequently, the velocity information needs to be provided somewhere. There are three possible approaches, each with different tradeoffs.

/**
 * Transforms a stationary point from source CRS to target CRS.
 */
DirectPosition transform(DirectPosition p);

/**
 * Transforms a point moving at the given velocity relative to the source CRS.
 */
DirectPosition transform(DirectPosition p, Velocity v);

/**
 * Gets the derivative of this transform at a stationary point.
 */
Matrix derivative(DirectPosition p);

/**
 * Gets the derivative of this transform at a point moving at the given velocity.
 */
Matrix derivative(DirectPosition p, Velocity v);

8.7.  Compact WKT

The Well-Known Text (WKT) format can be verbose because of repetitions. For example, if a file provides many CRS definitions, the same DATUM element may be repeated inside many of GEODETICCRS elements. The expansion becomes bigger when the file provides COORDINATEOPERATION definitions, because the same source and target CRS can be repeated many times and each repetition consumes tens of lines. Consider for example a concatenated operation:

CONCATENATEDOPERATION["A to D",
  SOURCECRS[A[…]],
  TARGETCRS[D[…]],
  STEP[COORDINATEOPERATION[SOURCECRS[A[…]], TARGETCRS[B[…]], …]],
  STEP[COORDINATEOPERATION[SOURCECRS[B[…]], TARGETCRS[C[…]], …]],
  STEP[COORDINATEOPERATION[SOURCECRS[C[…]], TARGETCRS[D[…]], …]]]

Figure 26 — CRS repetitions in concatenated operations

Each CRS (A, B, C, and D in above example) is defined twice. If some of those CRSs are derived or projected CRSs, then in addition to a CRS repetition, the operation parameters inside the CRS definitions are also susceptible to repetition in the COORDINATEOPERATION body. This repetition problem does not occur in the definition of a single CRS, which is the main purpose of WKT format. But it occurs and can become quite acute when defining a coordinate operation, or when defining more than one CRS in the same file.

GML can avoid this problem by defining each element only once, then referring to previous definitions using xlink. An example applied to CoordinateOperation was shown in Clause 5.6.1. Links exist also in YAML, but the WKT format and JSON (when no extension is added) lack an equivalent mechanism.

For reducing verbosity, the WKT format could be extended with aliases. This is not exactly the linking mechanism described in above paragraph, but it achieves similar results for the purpose of compactness. The WKT format could be extended by giving special meanings to the = and $ characters when they appear outside quoted text. Those characters are not valid in Unicode identifiers and cannot appear outside quoted texts according to the current WKT specification, so this proposal would not create ambiguity. The extension would introduce the following syntax:

<Unicode identifier> = <WKT definition fragment>
!! Restriction: the same Unicode identifier cannot be assigned twice in the same file.

Figure 27 — Definition of an alias

The WKT fragment can be anything provided that brackets or parenthesis are balanced. The occurrence of last closing bracket or parenthesis identifies the end of the WKT fragment. Then, the $<Unicode identifier> would be allowed to appear anywhere outside quoted text in the same file. The effect would be as if the <WKT definition fragment> associated to the Unicode identifier was copied inline.

The following example creates DEG, METRE, WGS84, LATLON, and WORLD aliases. The alias names have no special meaning and can be anything. An alias can use previously defined aliases (e.g. DEG in this example). The WKT element which is not an alias (when there is no = sign) is the main definition, which is GEODCRS["WGS 84", …] in this example:

DEG    = ANGLEUNIT["degree", 0.0174532925199433]
METRE  = LENGTHUNIT["metre", 1]
WGS84  = DATUM["World Geodetic System 1984",
           ELLIPSOID["WGS 84", 6378137.0, 298.257223563, $METRE]],
           PRIMEM["Greenwich", 0.0, $DEG]

LATLON = CS[ellipsoidal, 2],
           AXIS["Latitude (B)", north, ORDER[1]],
           AXIS["Longitude (L)", east, ORDER[2]],
           $DEG

WORLD  = AREA["World."],
         BBOX[-90.00, -180.00, 90.00, 180.00]

GEODCRS["WGS 84", $WGS84, $LATLON,
  SCOPE["Horizontal component of 3D system (…snip…)."],
  $WORLD, ID["EPSG", 4326]]

Figure 28 — WKT alias examples

The WGS84, LATLON, and WORLD aliases are not useful when defining a single CRS as shown above, but become more useful if EPSG:4326 is followed by the definitions of many other CRSs in the same file. Aliases can also be used for full CRS definitions, which is very handy for the CONCATENATEDOPERATION case mentioned at the beginning of this section.

This extension has existed for 7 years in Apache SIS (Clause 9). That implementation shows that this extension proposal can work, but it does not mean that an OGC standard would need to be identical.

8.8.  User-defined transformation registry

The EPSG database contains not only the definitions of thousands of CRSs but also contains transformations between various pairs of CRSs. Those transformation definitions are used by late-binding implementations of map projection libraries. A consequence of this architecture is that if a CRS has the following characteristics:

• is defined outside the EPSG database;

• is not a ProjectedCRS, DerivedCRS, or `BoundCRS`or is one of those types not defined in the EPSG database; and

• the CRS definitions are the only information (no `COORDINATEOPERATION[…] supplied),

then the software may not be able to perform coordinate transformations. Software will need some mechanism for completing the EPSG database with a user-defined registry of transformations.

A simple approach could be to encourage data producers to write COORDINATEOPERATION elements in Well Known Text (WKT) or GML format in a single text file. A single file could contain an arbitrary number of elements. Using the compact WKT proposal (Clause 8.7), the WKT file should be reasonably small and efficient to parse. Data producers would publish that file at some stable URL on their web site. CRS definitions would contain a link to this transformation registry. Software could then load that registry when first needed.

One issue is to decide how to specify the URL to a geodetic registry where map projection libraries can download coordinate operations. A possible location using current ISO 19111 and ISO 19115 models could be at the following path:

CRS.name → Identifier.authority → Citation.onlineResource

Where the OnlineResource would have the following property values:

• linkage = URL to the WKT file containing coordinate operations;

• applicationProfile = (to be determined); and

• function = coordinateOperationRegistry (a new code list value added in OnLineFunction). Alternatively, the function could be the existing download standard value if there is another way to specify that the purpose is to provide coordinate operations in WKT format.

The above proposal can be done without changing the ISO 19111 model. However, it would require a new keyword in the WKT format.

8.9.  User-defined operation methods

The definitions of Coordinate Transformations are restricted to the set of Operation Methods known to the implementation software. Each operation method needs to be associated to code in a programming language such as Java or C/C++; there is currently no way to define an executable OperationMethod in an implementation independent way. The best that can be done is to standardize the names of methods and parameters, then describe their use in a human-readable document such as EPSG guidance notes or TestBed 18 D023 Engineering Report. This approach works well as long as the standardized methods are sufficient. However, if a user requires formulas which are not included in the set of predefined operation methods, current OGC standards provide no help.

Implementation independent (in a given programming language) operation methods are technically feasible. Some programming languages have a concept of Service Providers which allows users to augment some aspects (services) of a library. For example ,in the Java language, ISO 19111 operation methods are represented by instances of the standard org.opengis.referencing.operation.OperationMethod Java interface (defined by OGC GeoAPI 3.0.2). If a GeoAPI implementation chooses to accept user supplied OperationMethod instances, then that implementation would put the following line in its module-info.java file:

module com.mylib {
    requires org.opengis.geoapi;

    uses org.opengis.referencing.operation.OperationMethod;
}

Figure 29 — Java implementation accepting user supplied operation methods

If a user named Alice defines her own operation method in a class named people.alice.MyTrajectoryConverter, then she can plugin that operation in the com.mylib library with the following lines in her module-info.java file:

module people.alice {
    requires com.mylib;
    requires org.opengis.geoapi;

    provides org.opengis.referencing.operation.OperationMethod
        with people.alice.MyTrajectoryConverter;
}

Figure 30 — Custom operation method added to a Java library

However, for making the above code useful, the OperationMethod Java interface currently defined by GeoAPI 3.0.2 would need to be extended. A key goal of the GeoAPI development work is to translate existing OGC/ISO abstract models to Java interfaces. GeoAPI does not add classes or properties that are not defined by other OGC/ISO standards, except for better integration with the target programming language. The properties defined by ISO 19111 OperationMethod class provide metadata about the coordinate operation but provide nothing for materializing the method in executable code. For the use cases described in this Engineering Report, the following method would need to be added to GeoAPI OperationMethod interface:

/**
 * Creates a math transform from the specified group of parameter values.
 *
 * @param  factory     the factory to use if this constructor needs to create other math transforms.
 * @param  parameters  the parameter values that define the transform to create.
 * @return the math transform created from the given parameters.
 * @throws InvalidParameterNameException if the given parameter group contains an unknown parameter.
 * @throws ParameterNotFoundException if a required parameter was not found.
 * @throws InvalidParameterValueException if a parameter has an invalid value.
 * @throws FactoryException if the math transform cannot be created for some other reason
 *         (for example a required file was not found).
 */
MathTransform createMathTransform(MathTransformFactory factory, ParameterValueGroup parameters)
        throws InvalidParameterNameException, ParameterNotFoundException,
               InvalidParameterValueException, FactoryException;

Figure 31 — Materializing OperationMethod to executable code in GeoAPI

This is demonstrated in the org.opengis.testbed.T18D025 Java demo in the OGC GitHub repository using an interface specific to the Apache SIS project (Clause 9.1) as a workaround for the absence of the above method in GeoAPI interface.

An alternative approach would be to specify formulas in CQL2, which supports trigonometric functions, among others. CQL2 may be sufficient for relatively simple coordinate operations, but more advanced coordinate operations may still need to be encoded in a programming language. For example, reverse map projections often use iterative algorithms, which require loops with non-trivial stop conditions. Coordinate operations may also need to load files, for example the trajectory in Clause 5.7. Coordinate operations are often complex enough to require the definitions of private methods, use of exception handling, etc.

8.10.  Reverse operation method

A coordinate operation sometimes needs to use the reverse of a coordinate operation defined in a geodetic registry. For example, “Mercator” is defined in the EPSG database as an operation method from geographic CRS to projected CRS. However, there is no registered operation for doing the coordinate operation in the reverse direction. In some cases, the reverse operation can be executed by using the forward operation with parameter values of opposite sign. But this is not the case of every operation method.

The legacy OGC 01-009 specification resolved this ambiguity by introducing an INVERSE_MT keyword in WKT format. This keyword means that the coordinate operation to apply is the reverse of the coordinate operation contained inside the INVERSE_MT[…] element. There is no equivalent feature in ISO 19111 (abstract model), ISO 19162 (WKT), or in GML format.

ISO 19162 discusses this issue in the WKT representation of concatenated coordinate operations section. The ISO approach is based on conventions on the CRS declaration order. For example, if the target CRS of a ConcatenatedOperation is the source CRS of the last step, the last step must be executed in reverse direction. This approach works in the context of concatenated operations. For the reverse direction of a standalone operation, a workaround may be to wrap the operation in a concatenated operation having a single step.

The Testbed 18 participants invented a scenario where coordinates are transformed from Voyage CRS to observatory CRS. The scenario was defined in that direction rather than the opposite direction in order to avoid the difficulty of expressing the reverse of EPSG:9837 (Geographic/topocentric conversions) operation method. For allowing scenarios in both directions, one of the following actions should be considered:

• introduce something equivalent to OGC 01-009 INVERSE_MT keyword; or

• document in a more prominent place how implementations should infer the direction of an operation method in the general case (not only in concatenated operations).

8.11.  Summary of proposed extensions

Extensions that could be done in a separate standard:

• add a CelestialBody class (Clause 8.1);

• new Extent subclass for regions at a distance from a celestial body (Clause 8.2);

• add new axis direction codes for declination and right ascension;

• if desired (it may not be the case), new code list values for axis directions such as planetocentricX, planetocentricY, and planetocentricZ (Clause 8.3);

• add an InertialCRS class and its InertialDatum dependency (Clause 8.4);

• add a MinkowskiCS coordinate system (Clause 8.5);

• compact WKT for definitions of complex transformation chains (Clause 8.7);

• standard location in the CRS metadata where to specify the URL to a transformation registry (Clause 8.8);

• add trajectory content type to GGXF format (Clause 5.7.1); and

• add relativistTransformation content type to GGXF format (Clause 5.7.2).

Extensions requiring modifications to existing standards:

• add an association from ISO 19111 GeodeticReferenceFrame to CelestialBody (Clause 8.1);

• add an association from ISO 19115 GeodeticExtent to CelestialBody (Clause 8.1.1);

• remove the EngineeringCS union from ISO 19111 standard (Clause 8.5.1);

• new WKT keyword in ISO 19162 for the URL to a transformation registry (Clause 8.8);

• materialize user supplied operation methods as executable code (Clause 8.9).

• allow source and target CRS in <gml:Conversion> (Clause 5.6.1.1) — possibly a GML bug;

• upgrade GML schema from ISO 19111:2007 to ISO 19111:2019 (Clause 5.6.1) — maybe already done in ISO 19136:2020;

• allow some GGXF content types to have a single CRS attribute instead of sourceCRS, targetCRS and interpolationCRS (Clause 5.7.2); and

• add a CoordinateOperation.transform(DirectPosition, Velocity) method (Clause 8.6.2).

Potentially useful features derived from previous standards:

• coordinate operation derivative as a Jacobian matrix (Clause 7.1); and

• WKT keyword for reverse operations (Clause 8.10).

9.1.  Demo application

For demonstrating the effectiveness of GML documents shown in this Engineering Report, a demo Java application is provided in Testbed 18 GitHub repository. This application uses Apache SIS, which provides a map projection library capable to read and write GML documents. Apache SIS supports the advanced features used in the VoyagerToObservatory.xml and PoseToRoboticArm.xml files such as DerivedCRS, CompoundCRS, PassThroughOperation, and ParameterValueGroup. The demo application also illustrates the extension mechanism presented in Clause 8.9 by expanding the set of operation methods known to Apache SIS. The operation methods added by the demo application (necessary for parsing the GML files of Testbed 18) are:

• trajectory to Conventional frame;

• to circular orbit (Spherical domain);

• time-dependent longitude rotation (equatorial CS);

• geographic/spherical conversions; and

• pose irregular series.

For this demo, the above operation methods are implemented with simplistic formulas, sometimes a trivial identity operation. But the formula complexity has no impact on the interoperability discussion.

After the above operation methods are made available to Apache SIS using the standard java.util.ServiceLoader mechanism, the following lines of code are sufficient for parsing the GML document and transforming coordinates. In this code, XML.unmarshal and GeneralDirectPosition are specific to Apache SIS. But the results — CoordinateOperation and DirectPosition — are implementation neutral GeoAPI interfaces.

// Standard Java
import java.io.File;
import javax.xml.bind.JAXBException;

// Implementation neutral GeoAPI
import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.TransformException;

// Implementation specific
import org.apache.sis.xml.XML;
import org.apache.sis.geometry.GeneralDirectPosition;

public class Demo {
    public static void main(String[] args) throws JAXBException, TransformException {

        CoordinateOperation op = (CoordinateOperation) XML.unmarshal(new File("VoyagerToObservatory.xml"))
        DirectPosition  source = new GeneralDirectPosition(x, y, z, t);
        DirectPosition  target = op.getMathTransform().transform(source, null);

        System.out.println("Coordinates relative to Voyager: "     + source);
        System.out.println("Coordinates relative to observatory: " + target);
    }
}

Figure 32 — Demo application reading GML and transforming coordinates