Computational geometry on the surface of a sphere (or an ellipsoid of rotation) requires an efficient and purposeful coordinate system: a method to represent each distinct surface location by a distinct set of numbers - or, preferably, by a single number. There are different methods by which this correspondence is defined, called geographic (geodetic, cartographic, etc...) coordinate systems.

This text introduces such system, specifically well‑suited for calculations on a digital computer. It describes the coordinate format and the transformations required for interactions with conventional latitude/longitude (φ, λ) based spherical or ellipsoidal coordinate systems. It is implemented in computer code of open‑source Nemo Library of Round-world Computational Primitives.

This text does not get into the details of significant advantages of replacing angular (φ, λ) spherical coordinates and trigonometry with vector (di, dj, dk) algebra as a method of choice for dealing with the planetary computational geometry when using digital computers. The topic is however addressed in the documentation of the Nemo Library and the commentary of API demo programs that can be found in the extensive collection of Nemo Library API examples. A simple "planetography" illustration of the continuity characteristics of the two coordinate systems might be a useful introduction to the topic.

Surface normal as the identifier of a spherical or ellipsoidal point

Conventional geographic coordinate systems use φ and λ for describing, by numbers, a position of a point (P) on the surface of a sphere or an ellipsoid. These two numbers represent the direction of a surface normal (n) in the point defined by its φ and λ coordinates. The same direction can be alternatively defined by its three normalized direction cosines (i, j and k) in the 3‑dimensional (x, y, z) Cartesian coordinate system, with the origin in the sphere's centre. (All as depicted above).

Since the sum of squares of normalized direction cosines is (by definition) equal to 1, in order to keep around a large volume of surface point coordinates, we need not record permanently all three of those: any one of them (when required to perform some vector calculation) can easily be restored by subtracting the sum of squares of the other two from 1, and taking the square root of the difference. It is worth noting that on a digital computer evaluating a square root is an operation many times faster than evaluating trigonometric (sine, cosine, etc.) functions. Numerical stability demands that the "dropped" and "restored" direction cosine be the one with the largest absolute value.

The range of normalized direction cosines is from -1 to +1. While the magnitude of any one direction cosine can be easily and quickly obtained from the other two, its sign is lost in the operation of squaring and square root computation. A method must therefore be established to additionally save two critical pieces of information: which of the three (i, j or k) direction cosine must be re‑calculated, and whether it is positive or negative.

This is achieved by recording another, quite succinct piece of information about the point: its domicile (digital) "continental plate".

Digital‑Continental ("DigiNental") Plates

If we consider the faces of a cube tightly encompassing the sphere (with a radius of 1.0) representing the planet, we can associate every point on the spherical surface with one of the six hexagon faces. The classification is fast and simple: if the point's normal first in the triplet of direction cosines (dc[i], di for short - the cosine of the angle between the normal and X coordinate axis) is the greatest (by the absolute value) and positive, the point is associated with the cube face that lies in X = 1 plane, if negative, with the face in X = -1 plane. This can be extended to points with other direction cosine maxima (positive or negative). We can further name the direction cosine to be dropped and re‑calculated as R, and the other two as P and Q, in the the cyclic order of x‑y‑z coordinate axes. For ease of visualization, hexagon faces can be named by their association with the capital geographical features of the planet, in a digital equivalence of the manner by which the geological "continental plates" are named. The term "digital continental plates" (abbreviated "digiNental plates" further in the text) will be used for this classification.

The summary of previous paragraph can be presented in a table:

And depicted on a stero‑cylindrical World map:

Retained direction cosine scaling

If the location of a point on the planetary surface (where the planet's radius is approximately 6.4 million meters) is defined by three normalized direction cosines in 64‑bit floating point register image in computer or disk memory, the spatial resolution of coordinates is in the order of about a nanometre; unnecessarily high for any application that manipulates large volumes of planetary surface location data.

As explained above, the normalized direction cosine with the greatest absolute value can be removed, to be re‑calculated when required. The greater of the remaining two can have no lesser value than -(√ 2 / 2) and no greater than +(√ 2 / 2)). A value in this range can be scaled down to the range that can be represented in an unsigned 32‑bit integer - i.e., between 0 and 2³²:

If two retained direction cosines are scaled down to an unsigned integer in the range between 0 and 2³², this will result in the decrease of spatial resolution to 5 millimetres on average, 15 mm in the worst case. Knowing the digiNental plate of the point provides three sine qua non pieces of information necessary to reconstruct "on the fly" 64-bit di, dj and dk normal vector:

1. Which of the three original direction cosines has been "dropped".
2. The sign (+ or -) of the direction cosine obtained by the √ operation.
3. Which of the two retained direction cosines (p and q) represent which of the original element (di, dj or dk) of the coordinate triplet.

The procedure is simple: first nibble (i.e., "half-byte") of the most significant byte (MST, see below under [Endianness]) is used to record the digiNental plate of the point. Following that, half‑bytes of two retained direction cosines are stored in a 64‑bit memory word - in most contemporary 64‑bit hardware architectures and programming languages represented as an "unsigned long integer". The last nibble is shared by the next‑to‑last two bits of two retained direction cosines. The very last two bits of each retained direction cosine are lost as "rounding error". This "nibble‑interlaced" packing schema results in an 8‑byte memory unit that can be called "UniSpherical 8‑byte point coordinate" (possibly abbreviated to "Us8 coordinate"). This method of "data‑packing" is depicted in the following diagram:

The "pseudo‑code" of the transformation of point coordinates from three direction cosines into Us8 coordinate format is as follows:

1. Determine the integer number of the point's digiNental plate by comparing the sizes and signs of its three direction cosines.
2. Store the plate number as the first nibble of the MSB of the Us8 coordinate.
3. As per the table of plates and R, P and Q direction cosines, scale P and Q into two 32‑bit unsigned integers.
4. Following the nibble used for plate number, store alternative nibbles of P and Q into the nibbles of 64‑bit Us8 word, in the decreasing significance of the two source and one target bytes. Use the last nibble for the two next‑to‑last pairs of P and Q bits. Ignore the very last two bits of P and Q.

In geodesy and cartography, the term "mapping" is commonly used for coordinate transformation from some geometrically more complex data domain to a simpler one, and "un‑mapping" for the inverse. Thus "mapping" is used here for the transformation of the 3‑D direction cosine triplet to the Us8 coordinate format, and "un‑mapping" for the inverse. In the case of UniSpherical coordinates, un‑mapping computations would be the inverse of the procedure outlined in the pseudo code given above, with the following minor addition.

Since, in the process of interlacing P and Q into a Us8 coordinate, the two least significant bits are truncated, a 1 (i.e., one bit in the least significant position) is added to both before scaling them back to -(√2 / 2) and +(√2 / 2) range - as this is the arithmetic mean of the truncated magnitude of P and Q vector components.

The above pseudocode is implemented in the open source Nemo Library in C language. Bitwise operations are implemented using "shifts"; while somewhat less "readable" than "C Bit-Fields", they tend to produce faster executing code and are likely easier transliterated to other programming languages.

The resulting coordinate (note the singular - a point location on the planetary surface in UniSpherical coordinate system is represented by a single number; not, as is the case for φ and λ, by two numbers) can be presented on an orthographic planetography animation:

In the above, each pixel's color shade (from light to dark) is commensurate with the numeric value of the planetary coordinate of the pixel center. The digiNental plate centers are marked and connected by a line that indicates their numeric order. The Plate‑to‑Plate shade discontinuity is least along the boundary of the plates that are next to each in their numeric order. The maximum discontinuity of that numeric order is, as expected, along the boundary between plate 1 (Africa) and plate 6 (Antarctica) - conveniently occurring along a great circle segment separating twp plates where South Atlantic meets the Antarctic Ocean.

UniSpherical 4‑byte coordinates

For most applications that deal with large volume of location data scaling down P and Q to 32 bits each achieves an almost ideal balance between the required coordinate resolution and the coordinate data volume. In addition, the sorting and searching in linearly arranged 64‑bit "atomic objects" is directly supported by existing hardware and software devices.

There exists however a class of location‑specific data objects (for instance, atmospheric condition measurements taken by ambulatory sensors) where data volumes can grow extremely large but the millimetric coordinate resolution is quite unnecessary.

Applications that deal with such data can add a parallel UniSpherical mapping and un‑mapping procedure (Us4), where P and Q are scaled to 16‑bit unsigned integers. The transformation schema depicted in the diagram above would change just a little: the middle two bytes in P and Q, and the middle 4 bytes in Us8 coordinate word would be removed.

With this, the coordinate requirement would be halved (from 8 to 4 bytes) and the resolution would decrease (from 5mm on average, to 15mm in a worst case for Us8, up to approximately 100 meters on average and up to 900 meters for a worst case for Us4).

In addition to significantly reducing the coordinate storage requirements, Us4 coordinate format provides a way to generate UniSpherical coordinate ordering examples on smaller scale maps - such as the following map of the Us4 linear coordinate order map centred on downtown San Francisco:

Note that in mapping as well as in un‑mapping, for both Us8 and Us4 coordinate variants, the order of bytes in 32 or 64 "word" (in most programming environments a 32 or 64-bit unsigned integer) is of significance when binary data is exchanged between different hardware environments. (See below, under "Endianness" for details).

Some final notes:

Endianness

If the UniSpherical coordinates are written to an external text file (typically using hexadecimal numbers with %016lx output print conversion specifier) the coordinates can be read on a system with different internal byte‑ordering ("endianness") with no adidtional consideration. However, when in the interest of minimizing the bandwidth use the coordinate data volume must be reduced as much as possible, and without the complication of introducing compatible compression protocols between sending and receiving system, it might be worthwhile to perform the transfer of coordinate data in binary form. In such instances the canonical endianness of UniSpherical coordinate 64-bit (8 byte) unsigned integer should be that of little‑endian; i.e., the task of reversing the byte order should be assigned to those system components running on (nowadays much less common) big‑endian hardware.

A Sphere or an Ellipsoid?

Normal vector components (or φ, λ angles) that are the source for transformation to UniSpherical coordinates can represent the surface point location on either a sphere or an ellipsoid. Best practice of geographical information processing requires that the system design includes a rigorous correspondence between the two, instead, as is often the case, of using spherical geometry productions on ellipsoid coordinates (for instance, determining the distance between two points using spherical great circle trigonometric calculation on ellipsoid point φ and λ).

Nemo Library provides a very fast, "near-conformal" sphere to ellipsoid mapping and un‑mapping transformations that use neither transcendentals nor series expansions. While there is nothing to prevent the system architect to transform geographical (i.e., ellipsoid) coordinates directly into UniSperical ones, it is advantageous include this step into the φ, λ to/from Us8 "computational pipeline". For details, please consult extensive documentation of the near‑conformal sphere (NCS) in the Library reference material.

Drift

After the initial truncation that occurs when the point coordinates are transformed from direction cosines triplet to Us8 format, further back and forth transformations, no matter how numerous, introduce no "location drift".

Execution speed of coordinate transformations

The following results were obtained by running gcc-compiled, non‑optimized C language code on hardware with BYTE UNIX Benchmarks (Version 5.1.3) characteristics:

Machine: x86_64 (x86_64) 12-core AMD Ryzen 5
12 CPUs in system; running 1 parallel copy of tests:
System Benchmarks Index Values        BASELINE       RESULT    INDEX
Dhrystone 2 using register variables  116700.0   57277073.9   4908.1
Double-Precision Whetstone                55.0       8381.5   1523.9

NCS to Us8: 250 million points in 7.554 seconds or 33.095 million points/second
Us8 to NCS: 250 million points in 5.945 seconds or 43.734 million points/second