The making of a Global Grid - Remembering my escape from flatland

Quaternary Triangular Mesh (QTM) is a spherical subdivision scheme for encoding vector geodata across a planet as recursive triangular subdivisions of an octahedron. This model of location and its hierarchical coordinate system evolved from experiments with a hierarchical raster data structure for encoding terrain relief named DEPTH (Delta Encoded Polynomial Terrain Hierarchy). This 2D pyramid- type structure encoded attributes (surface elevations) explicitly and locations (grid cell indices) implicitly. The paper describes how DEPTH evolved into QTM through a global discrete data grid called Geodesic Elevation Model (GEM), which more resembled QTM than a raster model. It used DEPTH to encode surface elevations in a forest of geodesic triangular quadtrees instead of in a planar rectangular quadtree. All three models were designed to make limitations of data accuracy and scale explicit. DEPTH and GEM capture elevations as ranges of values that decrease as the structure densifies, and perform limited interpolation. QTM captures vector data by encoding spherical 2D locations to the extent that their positional accuracy, certainty, or precision warrant, but did not use DEPTH. This paper is a memoir that summarizes the thinking that went into these data models and explores how properties and deficiencies of one led to the other. It does not present any breakthroughs or new applications. Rather, it documents inventions that influenced subsequent developments of discrete global grids, and might still do so.


Introduction
This paper is written as a memoir, recounting the explorations of the author who developed several variable-resolution data structures to store, evaluate, integrate, and retrieve geospatial data, some of which may have influenced how discrete global grids evolved and is written as a memoir.
The first two projects happened at the Laboratory for Computer Graphics and Spatial Analysis, a well-known skunk works in the Harvard Graduate School of Design (1965Design ( to 1992)).I joined it as a master's student in 1968 and left as a Research Associate in 1984.Subsequent developments took place while working in private industry and as a PhD student in Europe.I estimate that over 45 years I have devoted more than five man-years hatching, refining, coding, documenting, and promoting one of these schemes or another.Not only did my efforts not enrich me, a large corporation patented some of my key ideas knowing full well that they were mine.That can happen, but when you reach a certain age, you start to want others to know what you did with your life.Don't expect a global summing up.This is just a personal reflection.I feel unqualified to assess the impacts this research may have had on the discrete global grid (DGG) or other geospatial research communities, and so I shall not try to limn intellectual lineage.
My research involved creating, analyzing and displaying topographic and statistical surfaces.At the time, surfaces were almost always modeled as regular 2-dimensional grids of either integers or floating-point numbers.Such grids were often created by interpolating sparse data sampled across a geographic region or an abstract domain.For example, digitized spot elevations or contour lines on topographic maps could be gridded by interpolating heights in some fashion.
Demographic data could also be treated this way.In 1978, I generated a time series of population density maps for the US by interpolating population counts from all the decennial US Censuses at county centroids to a set of data grids in Albers projection.The grids were mostly about 330 columns by 195 rows, but at least a third of their cells were outside the conterminous US.They didn't receive any data, but they needed to be there because grids are rectangular.Such a waste of computing resources, I thought at some point, after realizing that empty cells needed to be stored, read in, occupy memory, and be visited by any program that manipulates that data grid.
In the late 1970s, online and offline storage was still pretty expensive.Big Data was measured in 100K units.Parsimonious programming practices motivated me to put data grids on a diet..I knew that gridded surfaces contain considerable redundancy and that their degree of autocorrelation was quantifiable and could be exploited to store elevations in a lot fewer than the usual 32 bits per value.
It occurred to me to sort and rank grid cell values in a linear array from which duplicate values were then removed to use as a lookup table.That would let me replace the data values in the grid with small integers that index into the sorted elevation values.
How small could that lookup table be?I knew that 25% of cells in my grids lay outside of the conterminous US, and so would be represented by one single value in the lookup table.The remaining cells might contain many redundant values as well, depending on the topography.But then, in a 2D array of floating-point numbers representing a continuous surface, what constitutes a redundant value?It seemed to be convolved with the precision of the input data.If the data source for surface heights has, say, 0.5m vertical precision, when gridded, there are likely to be more duplicated values than were they captured to the nearest 0.1m of elevation.Likewise, a grid with 1m-wide cells should contain more duplicated values than a 10m grid having the same vertical precision.So clearly, vertical and horizontal precision (and accuracy) are linked.
Around then, image pyramids [22] started being touted, mostly for sending gray scale images through slow communication channels (remember dialup?).Could I treat surfaces like images and encode them hierarchically?That would mean computing one elevation to characterize the entire surface, breaking that down into four quadrants, then into 16, and so on, until the entire grid is represented.Given that an m x m grid could contain at most m 2 values, and that an m/2 x m/2 grid would have m 2 /4 values, why store those values at any greater precision than needed to enumerate them?After all, they are not source data, only local averages.But then, when you store all those resolutions, you end up with 1.333… times as many values as you started with-hardly a savings.
But wait-how many bits do each of those levels need?Just enough to capture how quadrant elevations vary at each stage, I reasoned.The higher the resolution, the more bits would be needed.So Kelly Chan (now at Esri) and I proceeded to preprocess grids to find mean, median, and range of values, and then pick a "typical" value for the grid as a whole (e.g., its mean, median, or range midpoint).To decide when sub-quadrant values differed from the current estimate, we used a criterion or "step size" of 1/4 of the data range.Based on the difference between a quadrant's typified value and its parent's value, it would be coded as moving up, down, or staying the same.At each level, step size (like cell size) would decrease by half.These three possibilities could be encoded in two bits per quadrant, including a value to flag cells having no data.Total storage was 1.3333 x 2 bits x the number of grid cells, or 2.6667 bits per cell overall, yielding a compression factor of 12 for a single-precision grid.We named this encoding scheme DEPTH (for Delta Encoded Polynomial Terrain Hierarchy).We set out to program a proof-of-concept in FORTRAN.He provided bit manipulation routines and I built a data structure for the sequence of sub-grids that stored them in a linear array accessed as a heap.We tested the algorithm on two data sets, a rugged area called Dent de Morcles in the Swiss Rhone Valley (created by Kurt Brassel) and the Northern Gulf Islands near Vancouver Island (created by Tom Poker).We characterized and compared input and output grids and mapped them in 2D and 3D. Figure 1 shows results for four levels of detail of the processed Dent de Morcles grid from [18].
Our "delta-Z" compression strategy worked, but came at a cost: aliased values.They were less smooth than the source grids, and had "cliffs" in areas of high relief and plateaus that failed to capture details in areas of monotonous relief.We found that dividing the current step size by a factor less than 2.0 could avoid some of these issues (1.618, the Fibonacci ratio, worked pretty well) but other artifacts-such as not reaching local maxima or minima-could result.But by smoothing the output grid, we got much better fidelity [18] Why dwell on an old project having no relationship to global grids?It's because soon after I stopped working on efficient grid storage methods, I noticed that DEPTH could be adapted to model terrain for an entire planet across at multiple spatial resolutions, which became my next major project.

A Diamond in the Rough
An intriguing aspect of DEPTH was its potential to integrate data observations into a data grid, essentially by throwing darts at it and allowing them to come to rest at a specified level of resolution.Conventional data grids lack such variability, even if the accuracy of input measurements is known.Both cell size and the number of bits per cell are fixed in a data grid.DEPTH hierarchies provide finer-resolution cells if the accuracy of input data warrants, but need a different data structure.
Although I never processed much data with DEPTH, that possibility percolated, melding in my mind with visions of vectors and clouds of points.Might there be a way, I wondered, to store dense or scattered or gridded spot elevations at appropriate resolutions in a way that could fill in the gaps between the initial observations?And could such a data store represent locations anywhere on earth, not just in arbitrary local rectangles?This was a bit of a pipe dream, but elevation data for the entire world was already coming available in crude forms, such as coarse latitude-longitude grids or isolines digitized from small-scale maps.
Other geospatial researchers had thought along these lines around the same time.The first approach to a global DEM I heard of modeled the globe as a cube, dividing its six faces in a four-fold or ninefold way.Tobler and Chen [23] proposed this solution, which has endured in the literature.In 1983, however, the only data structure I knew of that could contain a global hierarchical terrain database was the octree, an eight-fold decomposition of a cube designed to encode volumetric data (as "voxels").Octrees seemed to me like overkill for encoding the spherical relief.Later, "sparse octrees" were developed that might have worked better.Today they find application mostly in games and flight simulations.At the time I felt that dealing with all the octree pointers to interpolate a surface stored as an octree would be too messy a proposition.
I no longer remember what made me consider a vector-based solution.One motivation may have been my irritation at data grids that represented the world in equirectangular (plate carrée) projections, which greatly distort the figure of the earth.Why not grid the globe in some way that preserves its form and topology and avoids awkwardness at and near its poles?For quite some time I had known about Buckminster Fuller's geodesic domes and also his concept called Vector Equilibrium (a solid object that breathes), a working model of which is colloquially called a "Jitterbug" [20].In fact, I owned several working Jitterbug models created by my friend Dennis Dreher, a geometric sculptor.This object has eight triangular faces hinged in pairs at their vertices in a way that lets them swing around.In its "closed" configuration it is an octahedron.When twisted, it smoothly opens up along its edges to form a cuboctahedron with six empty square faces.As it continues to rotate, the "bug" closes up again to form an octahedron whose faces align differently.A Jitterbug can oscillate between its "octa" and "cubocta" states with almost no force applied to it.
Three important properties of Jitterbugs to know: (1) The cube and the octahedron are dual polyhedra; each face of one corresponds to a vertex of the other; (2) In its cuboctahedral state, a jitterbug occupies more volume than in its octahedral state; (3) Connecting each cube vertex to the neighboring three octa vertices yields a 12-sided polyhedron (a rhombic dodecahedron) with diamondshape faces.
This shape caused me to consider modeling the earth as a Jitterbug, both because bugs assume a continuum of configurations having varying radii, and because its initial faces could be subdivided recursively into triangles.And so, goodbye Flatworld.
Here's how a Jitterbug could host a global terrain model: (1) Imagine a wireframe cube with 6 green vertices.Elevations must be specified for the 14 GEM vertices.The four equatorial vertices of the octahedron could be set to the equatorial radius of the earth (~6378 km).Its two polar vertices elevations would likewise be set to the polar radius of the earth (~6356 km).Elevations for the eight cube vertices could be set at the mean volumetric earth radius (~6371 km).As this set of elevations defines an earth ellipsoid, any ellipsoid model could be used to specify them.What happens when you toss a 3D geographic location into GEM?Essentially, it traverses down to occupy a vertex of small triangle that befits its horizontal and/or vertical accuracy and is gobbled up.Its 3D coordinates can be tossed away because they can be reconstructed by an inverse traversal when retrieving the point.
Users of a GEM database would typically retrieve topographic data for a limited region.Besides its extents, they could also specify the accuracy they want the data to have or the map scale at which it should be displayed.GEM would oblige with sets of points or triangles.Retrieved points could be subsequently triangulated, if desired.Users could be confident of the accuracy of the terrain data returned to them, assuming the data's accuracy was respected when it was encoded into GEM.
When encoding a geographic point (presumably given as a 5-tuple: latitude, longitude, elevation, horizontal accuracy, vertical accuracy), it would enter one of the octants and be assigned its code (1)(2)(3)(4)(5)(6)(7)(8).Next, the closest octant vertex to the point is triangulated with the two cube vertices the point lies closest to, forming a triangle that connects the cube and octahedron.The next triangulation identifies which of nine sub-triangles of the octa face the point occupies, and its index (1-9) is appended to the GEM code.This procedure, called trilocation [16,17], recurses to progressively fix the point in the triangular hierarchy.
Triangle identifiers are only computed at odd-numbered levels, that is, for subdivisions of octahedral faces.At every step, sub-triangles are numbered 1 through 9, in an arbitrary order.In one of several possible GEM numbering schemes, the area around Wall Street and Ground Zero in lower Manhattan is covered by a triangle having the address 698258464, representing eight subdivisions of octa face 6 and which occupies about 1.5 km 2 .

Keep the Structure Simple, Stupid
After publishing my only paper about GEM (based on my AutoCarto Six 1983 presentation), I busied myself otherwise and then soon left academia.GEM was incarcerated in folders in my bookcase for five years and four domiciles before I revisited it, when working at a computer company that was entering the GIS marketplace.Wouldn't it be great, I thought, if GEM could be built into the new GIS?
Dusting off my paper and notes, I considered how to turn this odd duck of a data structure into a GIS game-changer, as we liked to say.Now more sensitive to public perception than I had been at the lab, I had to admit that the terms like "triacon breakdown of a cuboctahed-ron" (aka "vector equilibrium") might be a tad off-putting to paying customers.And as some details of GEM confused even me, I decided to recast the model to one I could communicate in a more user-friendly fashion.I did, but my work never made it into the product for complicated reasons.
My first thought was how to make the cuboctahedron go away.I sensed that little had to change to morph the model into a standalone octahedron.I also realized that hierarchical tessellation could help manage all sorts of geodata, not just terrain.After all, GEM does encode geographic locations globally in a scale-specific way, and that capability could serve a variety of purposes.In that light, topographic elevation is just one attribute of location, even if cleverly encoded.
As GEM codes only represented partitions of the octahedron, not the cube, maybe dual tessellations weren't needed.It would be natural to recast GEM as a set of eight triangular quadtrees having triangular facets.In Bucky-speak, this arrangement is called an alternate partition of a triangle, as opposed to GEM's triacon partitioning model.There are others that I ignored.Divided Spheres, by Edward Popko [20] tells all you ever need to know about partitioning spheres into polyhedra.
The name I settled on for the new model was octahedral quaternary triangular mesh (O-QTM), which described the geometry of the breakdown of the octahedron, a bit of a pun on UTM [17,18].
A QTM grid develops from an octahedron by bisecting octant edges at their midpoints and connecting those points with new edges to form four triangles.Figure 3 shows the first two levels.Only the octant edges follow great circles.All subsequent edges that partition facets follow small circles.These three orthogonal sets of edges, one set of which follows parallels of latitude, define the tessellation.Thus the shapes and sizes of the resulting facets are unequal.I worried about that for a  QTM facets occupy a forest of eight triangular quadtrees.Trilocating in QTM does not converge on locations as quickly as for GEM facets because facet size decreases by one-quarter versus one-ninth.Also, instead of the strings of decimal digits that are GEM location codes, QTM uses strings of two-bit quaternary numbers (0-3), called QTMIDs.For example, the location of my office at the University of Zürich, measured from a Swiss topographic map, was 47˚ 23' 48" N, 8˚ 33' 4" E. It has the level 18 QTMID 1133013130312301002, which is accurate to within 60m.The initial 1 is the octant code, and requires four bits; the rest use two bits.It geocodes a triangular area of about 900 m 2 in 40 bits.In a conventional geospatial database, storing a point location at a specific level of precision would take three 32-bit or 64-bit words; a longitude (or x) value, a latitude (or map y) value, and the measurement's positional accuracy (given as a spherical or planar distance).
A 32-bit QTMID can distinguish locations within a kilometer of each other, which is adequate for mapping at 1:1 million and up.Storing QTMIDs in 64-bit words makes possible 30 levels of detail.Such a QTM grid has 1cm resolution, which is suitable for mapping at 1:20 scale.
Trilocating a geographic coordinate point to obtain a QTMID for it requires traversing one of eight quadtrees.Given the resident octant, the process appends QTM codes to its ID using an efficient geometric algorithm called TRILOC, which computes spatial relations in Manhattan space.That space is a world map projection called Zenithial OrthoTriangular, or ZOT for short.ZOT projects the eight faces of an octahedron to eight right triangles occupying a square.Figure 2b shows a cubocatahedron in ZOT space.Within it, the octahedron's vertices are shown by red diamonds.Computations in ZOT space are cheap because no squares or square roots are needed to calculate distances.The TRILOC algorithm recursively converges on a QTM facet, stopping when the facet's edge length would be smaller than the error or accuracy tolerance specified for the point being encoded.The algorithm saves its stack to avoid repeating calculations when the next point lies near to the one just encoded, which is typical.The TRILOC algorithm is documented in [7] and [4].
As it traverses the hierarchy, TRILOC enumerates QTMIDs, one digit at a time.As there are 16 possible ways to enumerate four things, I could have arbitrarily picked one of these permutations.However, as orientation of QTM facets changes (the central ones flip up and down), I realized they should have a more coherent numbering system that expresses their spatial relationships.The one I chose codes QTMIDs such that all point-neighbor facets end with the same digit (1,2, or 3), and the IDs of all central facets end with 0. Figure 4 is an equicylindric map of the first three QTM levels (512 of them across the planet) color-codes facets according to the last digit in their QTMIDs.
This map reveals how QTMIDs cluster around vertices to form non-overlapping hexagonal regions (except at the six octa vertices, where they are diamonds).The triangular areas (yellow) are central facets that are always coded 0. I named the hexagonal clusters attractors (a chaos theory term) because their central nodes attract IDs in their vicinity [17].QTMIDs systematically reflect the basis numbers of the initial octa vertices, the digits 1, 2, and 3, using 0 to label central facets down the line.These numbers identify axes of the octahedron.Its North and South poles are numbered 1.The vertex located at 0° N, 0° E, and its antipode are numbered 2. The other two equatorial vertices are numbered 3. Vertices that define octa edges thus have different basis numbers, either 1-2, 1-3, or 2-3.Bifurcating an edge generates a new vertex at its midpoint, which receives whatever basis number edge's endpoints don't have (3, 2, 1, respectively).The formula for assigning a basis number (BN) is: BN new = 6 -(BN from + BN to ).The new node becomes an attractor associated with that basis number.
Attractors (whose IDs systematically alias the QTMID of one of the facets constituting them) proved useful for my application, generalizing vector map data across a range of map scales [3].This became the focus of my doctoral graduate research at the University of Zürich (1995-1998), after letting QTM lie fallow for several years [4].There I re-implemented code for QTM spatial transformations in the C language, and developed algorithms for generalizing vector map objects represented as lists of QTMIDs.

Getting to the Point
Map generalization was not the main motive for developing either GEM or QTM, although I realized that scaling map representations could indeed be a useful application of hierarchical coordinates.My primary aim was to identify locations anywhere on a spherical body compactly and in a way that expresses their locational accuracy or certainty.
In 2000, I took part in the first International Conference on Discrete Global Grids, and was thinking about QTM as a geocoding convention that would enable efficient and reliable validation, integration, and exchange of geographic data.The name I gave for this larger concept was a pun, Loc8.The term Loc8 identifier took the place of QTMID [1].Loc8 identifiers can build in qualitative metadata for geographic features at the coordinate level, allowing each point primitive to express its uniqueness.
If features from different sources or layers coincide, coding their coordinates as Loc8 enables them to be aliased at some level of detail.If the Loc8 codes do not match as precisely as they should, then perhaps their original representations, scales, or methods of capture differed.In some cases, Loc8 can assist merging different representations of topological primitives (points, arcs, polygons).
Loc8 codes are QTMIDs plus fields of metadata.Stored as 64-bit binary words, QTMIDs of 30 levels or less leave bits available to code additional metadata.The schema in Figure 5, from [4,8], shows 14 high-order bits containing seven 2- bit codes qualifying the QTM coordinate in the low-order 40 bits (level 18, 40m resolution).Six zero bits delimit the seven qualifiers from the QTMID.Two more QTM levels or qualifiers could occupy that null field if desired.This example shows point metadata for use in map generalization.Each qualifier is a trinary number coded as 2-bit fields (01, 10, or 11; 00 serves as a delimiter, not a legal qualifier).The table indicates how these qualifiers hint at roles a point and features that use it may play in a cartographic context.Clearly, qualifiers with entirely different meanings could be substituted, and that they need not be two bits as long as they are nonzero.Essentially, Loc8 encodes ontologies at the coordinate level.The bit fields can be made wider to accommodate more than three codes.
If something like Loc8 adopted as an international geospatial data standard, it could assist data exchange and integration while providing inherent certification of positional accuracy.One way in which Loc8 could promote data exchange and reuse is shown in Figure 6.
How geospatial databases represent geometry, topology, and define features and metadata can vary greatly.When geospatial data are transferred to different database formats, incorrect, lossy, and tedious data integration can result.Data exchange formats such as DIGEST and SDTS exist to facilitate integration, but success requires considerable external metadata.The schema shown here only describes exchanging coordinate data and metadata for it.Other data content (such as simple linear and areal objects) would require its own metadata to achieve interoperability.More complex objects can be built from these primitives.One could code other compact, hierarchical location codes from DGGs other than QTM using Loc8 identifiers.The grid could define triangles, rectangles or hexagons.They only need facets that completely cover the globe and can be regularly decomposed into similar figures that nest (properly or not).However, when a DGG defines facets having different tessellations, shapes, orientations, or hierarchies than some other system uses, transformation functions would be needed to be able to exchange data between models.Converting DGG location codes to spherical coordinates and from them to the other system's location codes may be required if more direct transfer functions cannot be devised.
There are many types of DGGs, and no two of them code locations the same.This may be inevitable, but it is also quite inconvenient.Given that no particular DGG can serve all geospatial applications, methods to convert between them need to be developed and made available, perhaps by a standards organization such as OGC (probably through its Discrete Global Grid Systems working group) or ISO (through ISO/IEC 18026:2009, Information Technology: Spatial Reference Model).Loc8 (which is not patented), perhaps served from middleware, could be used as an intermediate representation as well as a format for data exchange, somewhat along the lines of aspects of OGC's netCDF (raster) or Simple Feature Access (vector).

Conclusions
A longstanding interest in statistical and terrain surfaces and the need to represent them efficiently launched a series of investigations beginning in 1979 that led to two early digital global grid models.This paper has explored the paths taken to create this body of unfinished work, which for all its deficiencies inspired subsequent developments in that new field.Both the lessons learned and the current state of the art underscore ongoing needs to standardize hierarchical geospatial referencing models and methods.If incorporated into geographic information systems, hierarchical grids can do much to document positional data quality and facilitate data exchange.Advances in this area can foster new spatial analysis techniques and cartographic representations, and help to make global grids serve more diverse applications and eventually become more globally interoperable.

Figure 1 :
Figure 1: Four stages of DEPTH elevation encoding for the Dent de Morcles region of France (2) Surround it with a wireframe octahedron with 8 red vertices.(3) Place the vertices of the cube and the octahedron at the same distance from the origin.(4) Rotate the cube so that it's green vertices are in the center of the red triangles.(5) Connect each green vertex to each of the three red vertices surrounding it.This object, shown in Figure2, defines two sets having 24 triangles each: (a) 8 sets of 3 triangles having one green and two red vertices, and (b) 12 sets of 2 triangles having two green vertices and one red vertex.One could fully define the figure (a rhombic dodecahedron) with either set of triangles by eliminating either the red-red edges or the green-green edges.

Figure 2 :
Figure 2: The dual tessellations of the Geodesic Elevation Model (GEM) I chose the second set of triangles (having 2 green vertices and one red vertex) as the basis for what I called geodesic elevation model (GEM).I preferred it because those triangles are more nearly equilateral.By expanding or contracting the size of the octahedron or the cube, the initial dodecahedron would flex along the edges of the cube or the octahedron, respectively.Elevations must be specified for the 14 GEM vertices.The four equatorial vertices of the octahedron could be set to the equatorial radius of the earth (~6378 km).Its two polar vertices elevations would likewise be set to the polar radius of the earth (~6356 km).Elevations for the eight cube vertices could be set at the mean volumetric earth radius (~6371 km).As this set of elevations defines an earth ellipsoid, any ellipsoid model could be used to specify them.