WO2018136996A1

WO2018136996A1 - Digital earth framework and indexing schema to extend any "surface" discrete global grid system to construct interoperable "volumetric" discrete global grid systems

Info

Publication number: WO2018136996A1
Application number: PCT/AU2018/050042
Authority: WO
Inventors: Matthew Purss
Original assignee: Pangaea Innovations Pty Ltd
Priority date: 2017-01-24
Filing date: 2018-01-24
Publication date: 2018-08-02
Also published as: AU2018213401B2; AU2018213401A1

Abstract

A digital Earth framework and associated indexing schema that enables a given two-dimensional polyhedral representation of the Earth's surface, (a "surface" Discrete Global Grid System (DGGS)), to be extended to construct a three-dimensional "volumetric" DGGS. The invention comprises extending the cells of a given "surface" DGGS tessellation in Earth Centred Earth Fixed (ECEF) coordinates from the reference ellipsoid (e.g. WGS84) both inward to the centre of mass of the Earth and outwardly to an arbitrary outer ellipsoidal shell defined above the reference ellipsoid. The radial dimension is divided into equal-volume ellipsoidal shells that are each tessellated using a chosen "surface" DGGS. The equal-volume voxels of the "volumetric" DGGS are constructed by joining the vertices of each equal-area "surface" DGGS cell defined on successive ellipsoidal shells. Each "volumetric" DGGS voxel is referenced at its centroid using a composite indexing schema.

Description

Title of Invention : DIGITAL EARTH FRAMEWORK AND INDEXING SCHEMA TO EXTEND ANY "SURFACE" DISCRETE GLOBAL GRID SYSTEM TO CONSTRUCT INTEROPERABLE "VOLUMETRIC" DISCRETE GLOBAL

GRID SYSTEMS

Technical Field

[0001 ] The embodiments described herein relate to systems and methods for storing and analysing volumetric geospatial data associated with three-dimensional geospatial locations. In particular, the methods and systems pertain to the extension of systems and methods for providing discrete global grid systems (DGGS) designed for storing data at the surface of the Earth to include the three- dimensional storage of volumetric geospatial data observed within the interior of the Earth, at the Earth's surface, within the Earth's atmosphere and into outer space referenced to the location of the Earth's centre of mass.

Background Art

Storing and Representing Data Using Discrete Grids

[0002] Spatial information is commonly represented on a digital device using a

tessellation of discrete cells to represent a regular division of space. These representations of data are often referred to as "raster images" where each discrete cell is known as a "pixel" and these pixels are arranged in a flat (or planar) two-dimensional rectangular array of rows and columns. The value of each pixel represents the average interpolated value (for the area covered by that pixel) of the set of spatial observations used to construct the raster image.

[0003] Planar, regular rectangular grids are commonly used because of their direct correlation and compatibility with planar Euclidean geometry and the associated simplified matrix and array operations that can be leveraged from the grid structure. It is important to note, however, that raster grids and images are simply model representations of an original set of spatial observations. [0004] Other model representations of spatial information include vector objects; which are interpretations of spatial features using geometric shapes (points, lines and polygons) that are derived mathematically from a set of vertices or nodes referenced by spatial coordinates. Rectangular Cartesian coordinate spaces are often used as reference frames to describe point locations and vector quantities associated with these geometric shapes, which enables a common reference frame for the representation of vector objects as raster images, and vice-versa. Although, choices need to be made in any such transformations regarding and reconciling differences in spatial resolutions between different raster and vector objects. These choices can lead to increased error and uncertainty in the information represented by a particular discrete grid.

[0005] In three-dimensions, it is commonly the case that the planar two-dimensional rectangular grid is merely extended into an orthogonal third dimension to form a regular three-dimensional Cartesian grid that represents the spatial information as a three-dimensional matrix (or array) of rectangular prisms; each representing the three-dimensional average of the interpolated values of the original set of spatial observations. This extension also provides for a simple dimensional extension of the same matrix and array operations that are applied to a two- dimensional rectangular grid to also be applied in the three-dimensional rectangular grid.

[0006] There are numerous assumptions in this approach, particularly regarding the spatial reference frame of the original observations, and it is commonly assumed that the spatial coordinates of the set of observations on the surface of the Earth lie on a "flattened" Earth Surface - allowing a one-to-one translation into the three-dimensional rectangular grid space to both represent data and apply computational processes to the data. At local scales, and in the absence of variable topographic relief (e.g. large hills and valleys), this assumption is usually reasonable, given the inherent uncertainties in the measurement of the original spatial data observations. At larger lateral scales; however, this assumption breaks down due to the effect of the curvature of the Earth's surface. In these cases, other grid representations of the spatial information are required to more accurately represent the spatial relationships between a given set of spatial observations. Failure to correctly account for the curvature of the Earth for spatial analyses over large regions can lead to significant and costly errors in model outputs.

Polyhedral Tessellations of The Earth

[0007] There are other methods of dividing a given space to represent spatial

information as a discrete grid. One of these methods includes the use of convex polyhedrons (e.g. Platonic Solids) to tessellate the surface model of the Earth and construct a Discrete Global Grid System (DGGS). A "surface" DGGS is a spatial reference system constructed by mapping the faces of a platonic solid (i.e. tetrahedron, cube, octahedron, dodecahedron or icosahedron) to the surface model of the Earth to form a tessellation of cells representing equal areas on that surface. This initial tessellation is then refined by dividing each cell by a fixed areal refinement ratio to form a corresponding number of equal area child cells. This refinement process can continue infinitely to form infinitesimally small equal area cells (although for most current use cases will not require cell refinement smaller than ~1 cm²).

[0008] DGGS, unlike commonly used rasterization methods, are not based on a fixed planar rectangular grid representation of the Earth. The tessellation directly represents the curved surface of the Earth. The cell shape is also not fixed, but can be triangular, square, hexagonal or a mix of these shapes; usually, related to the chosen polyhedron used to construct the "surface" DGGS. This enables different styles of "surface" DGGS to be constructed to suit particular use cases, and/or data representation needs.

[0009] Most importantly "surface" DGGS provide not only an equal-area tessellation of the Earth's surface, but also a fractal-like hierarchy of discrete global grids with resolutions ranging from infinitesimally small to global in scale. This provides a common framework that simplifies the integration of spatial information across multiple scales and resolutions. It also enables the direct integration of both raster and vector objects without requiring the resampling, or vectorisation, of either object. A "surface" DGGS treats all data stored within it merely as a set of spatial observations that are assigned to individual cells in its hierarchy of resolutions. Whether a piece of spatial data is part of a raster dataset or a vector dataset is merely contextual information that may be used with other contextual and metadata information during analytical processes. Furthermore, a "surface" DGGS facilitates the spatial encoding of observation and measurement uncertainty, where an observation may be stored in the "surface" DGGS at a resolution that represents the level of lateral uncertainty with the spatial location of that observation. This greatly simplifies the consideration of uncertainty in spatial operations.

[0010] For many use cases requiring the integration of spatial information at the

Earth's surface, a "surface" DGGS will suffice. However, there are many applications where three-, and four-dimensional integration of spatial data is desired, or even required (e.g. climate models). It is common for these types of applications to use a spherical model of the Earth's interior, its surface and its atmosphere, tessellated into three-dimensional pixels, often referred to as voxels, using octree and similar methods to subdivide the volume of the spherical Earth model into regular units. These octrees are commonly constructed using either cubes or octahedrons to tessellate the spherical space and often use hierarchical indexing schemas derived from the topology of each individual face of the base polyhedron - not too dissimilar to common indexing methods used by many "surface" DGGS. While spherical models of the Earth are a reasonable proxy for the shape of the Earth, and also simplify the process of producing equal-area polyhedral tessellations, there is a significant level of spatial distortion between spherical and ellipsoidal models of the Earth. The latter being a more precise representation of the Earth's surface. Use cases requiring high precision, or accurate analyses at multiple scales (from molecular to global) should use ellipsoidal rather than spherical "volumetric" DGGS implementations.

Addressing Cells/Voxels Using Space Filling Curves

[001 1 ] Space filling curves are a common method to efficiently reference and

organise spatial data on a computer (or digital device) such that the data is stored in close association with nearby data. When stored in the order of the locations along the space filling curve, nearby data will also occupy nearby physical memory locations. This results in highly efficient storage, search and retrieval of spatial information. [0012] There are numerous types of space filling curves that have been published in the literature, each with their strengths and weaknesses. Of these the Hilbert and Morton (or Z) space filling curves are examples that are commonly used for referencing spatial information; Morton curves are often favoured over Hilbert curves because the mathematical formula that defines the path of the curve is less computationally expensive than that for the Hilbert curve; and thus, Morton curves are often more efficient in practice. Prior art of the application of two- and three-dimensional space filling curves is extensive; however, most applications use regular planar Cartesian realisations. Most commonly, the ordinate axes of a model space (or discrete grid) are converted, or mapped, into an integer set of subdivisions which often form the numeric ordinate references that are used to construct the space filling curve indices. There are some prior art examples of latitude and longitude angular coordinates being used in bitwise operations to encode a set of spatial (latitude/longitude) coordinates into an integer space filling curve reference. This invention leverages and extends this approach to encode the three-dimensional geodetic location of a spatial observation as an ellipsoidal coordinate space (i.e. latitude, longitude, radius) realisation of a Morton space filling curve.

[0013] Other prior art methods of encoding discrete cells and voxels along a space filling curve involves the hierarchical ordering of cells/voxels used to subdivide a discrete space and using individual space filling curves to reference the cells of each successive level of resolution.

Addressing Cells/Voxels Using Hierarchical Indices

[0014] Hierarchical indices are used extensively in prior art that underpins the

development of multi-dimensional relational and spatial databases to enable efficient data search and retrieval. In many cases these spatial data

infrastructures use a division of the Earth's surface achieved directly through the use of equal subdivision (or graticulation) of latitude and longitude angles - producing unit cells that are not equal in area or shape and that have no direct hierarchical topology.

[0015] The more sophisticated prior art spatial data infrastructures, such as DGGS, use balanced equal area subdivision methods (such as quadtrees and fractal subdivision) to produce equal area cells with an indexing schema that is inherently hierarchical and multi-resolutional. Extending to the three-dimensional case, the common approach of most prior art is to use octrees to construct the hierarchical topology and indexing schemas for voxels.

[0016] The levels of refinement of a "surface" DGGS also provides a form of

hierarchical reference that can be used as part of an indexing schema for cells and voxels of two- and three-dimensional grid infrastructures. Some prior art DGGS implementations (e.g. the rHealPIX DGGS) use the length of the cell index to reference the refinement resolution of cells, and thus indicate the position of a cell within the hierarchy of the DGGS. This invention combines a novel three- dimensional ellipsoidal space-filling curve reference with a reference to the resolution of a "surface" DGGS to construct a composite hierarchical indexing schema for voxels that is independent of the cell indexing schema of the

"surface" DGGS.

Systems That Use Spatial Data Representations

[0017] Spatial data are commonly represented on a computer using a Geographic Information System (GIS), which is a computer software application that enables the creation, storage, visualisation and analysis of georeferenced spatial data. A GIS may be configured to work with raster data, vector data or both, although, commonly a GIS will treat raster and vector data differently. A GIS which is based on a DGGS is able to work with both raster and vector data in the same way.

[0018] While some GIS are configured for two-dimensional spatial information (i.e. horizontal surface x/y or longitude/latitude spatial locations), others are also configured to work with three-dimensional information, either through contextual metadata means (e.g. spatial representation of horizontal coordinates coupled with a colour-scale representation of a given piece of information, such as height/elevation, associated with that location; these are 2.5 - dimensional systems) or full three-dimensional volumetric representation of spatial data.

[0019] The visual representation of multi-dimensional spatial information on a

computer using 2-, 2.5- and 3-dimensional GIS applications requires a translation of the geographic locations of spatial data to the two-dimensional digital coordinates of a computer monitor to enable both qualitative and quantitative analysis of data.

[0020] The value of spatial information is often in the additional contextual

information of an observation, not merely its spatial location. This contextual information may be descriptive information associated with a spatial location, such as place names or cadastral boundaries, or physical measurements of phenomena and processes of the Earth acquired from sensors and sensor networks.

Spatial, Digital and Geodetic Transformations

[0021 ] Like paper maps, a computer monitor is a flat two-dimensional medium. This necessitates the digital transformation of spatial information from geospatial coordinates to the two-dimensional digital reference frame of the monitor for visual representation of spatial information. In many GIS applications various geodetic transformations are applied to spatial information to produce planar map projected data first which is then transformed to the digital reference frame of the monitor for display. Often a similar set of geodetic transforms are applied for three-dimensional analysis systems with the vertical dimension represented by the map projected height/elevation information to produce model spaces that are planar Cartesian rectangular prisms. In order to integrate, analyse, model and interpret spatial information it is often necessary to apply geodetic

transformations and map projections to translate spatial data into a common reference frame. For many non-DGGS spatial reference systems a further step of resampling raster and vector data to the same resolution is also required to facilitate spatial analysis and data fusion.

[0022] The use of geodetic transformations and map projections simplifies the

process of applying digital transformations to represent both two- and three- dimensional spatial information on a computer (i.e. from one Cartesian reference frame to another); however, it can be problematic for precision spatial data integration, analysis and fusion. This is because many geodetic transformations and map projections result in some level of degradation in the accuracy and precision of spatial coordinates leading to increased spatial uncertainty and error in subsequent analyses. This uncertainty can accumulate with successive applications of geodetic transformations and map projections. The use of DGGS technologies minimises this issue by negating the requirement to apply map projections. DGGS based applications may require the an initial application of geodetic transformations on data import if the spatial data being imported was acquired using a different geodetic datum to that of the DGGS.

Summary of Invention

[0023] The object of this invention is to provide an efficient and flexible method

whereby any "surface" DGGS implementation may be extended to construct a three-dimensional "volumetric" DGGS implementation that encompasses the interior, surface and atmosphere of the Earth including outer space extending to an arbitrary maximum radius defined by an integer multiple of the Earth's volume and whereby each DGGS voxel is referenced by an indexing schema that is independent of the indexing schema employed by the "surface" DGGS used to construct it. This object is solved by a method according to claim 1 and a volumetric Digital Earth Information System (vDEIS) according to claim 6, claims 2 to 5 refer to specifically preferred realizations of the inventive method, claim 7 refers to the realisation of a distributed network of vDEIS made possible by the spatial data storage, search and retrieval mechanisms described by claim 6 that leverage the spatial indexing schema proposed in claim 5.

[0024] The invention realises an efficient method to construct three-dimensional volumetric model spaces that are scalable from infinitesimally small to global scales. The method of indexing provides a way to map three-dimensional geospatial data to a one-dimensional array representation on a computer data storage device that preserves the local spatial relationships of the data.

[0025] The invention provides for multiple distributed "volumetric" DGGS

implementations to be constructed from one, or a number of different styles of, "surface" DGGS while maintaining a common indexing schema that enables for efficient storage, search and transfer of geospatial data within and between each "volumetric" DGGS instance. Technical Problem

[0026] The Earth is a dynamic planet with a myriad of interrelated systems and

processes that range from microscopic to global scales. These systems impact on the human experience, just as we impact on them. In our thirst to understand these systems we are acquiring more and more data across an increasing range of scales. And to answer the increasing number of questions we ask as our understanding of the Earth and her systems increases we need to be able to convert this data into useful information; not just on the surface but within the interior of the Earth, and its atmosphere and surrounding space as well.

[0027] Issues such as climate change, food and energy security, epidemiology, natural resource and ecosystem management, migration and urban development - to name but a few - require the integration and analysis of geospatial data and information from multiple sources at multiple scales, in multiple dimensions and with an increasing number, and sophistication of, sensor networks.

[0028] As the volume, velocity and variety of spatial data increases the ability to

rapidly and efficiently produce meaningful information from the data we collect is becoming increasingly challenging. With the emergence of the 'Big Data Era' it is now clear that conventional GIS approaches, using layers of planar map projected data that require repeated rasterisation, vectorisation and resampling to facilitate spatial analyses, are not scalable to sufficiently answer the questions we are asking within the timeframes, data volumes and across the spatial scales required.

[0029] Integrating spatial data from multiple, often disparate, sources across multiple scales, resolutions and measurement precisions is a complex and non-trivial exercise. Conventional treatment of data in a standard GIS context produces raster and vector datasets that are internally consistent but require resampling in order to integrate them with additional spatial datasets. Often choices must be made that may degrade the spatial accuracy of one or a number of datasets in order to facilitate integrated spatial analyses. This repeated and often necessary transformation and resampling of data adds an additional processing cost to analyses that, while acceptable at small scales and volumes of data, becomes increasingly unworkable as the scales and volumes of data increases. Solution to Problem

[0030] In recent years the maturity of DGGS technologies along with high

performance and cloud compute environments has provided an alternative approach that solves many of these issues. DGGS provide a common analytical infrastructure that avoids the necessity to account for the distortions created by map projections of individual datasets; provides a fixed, equal-area grid framework that accurately represents the surface model of the Earth and consistently scales from infinitesimally small to globally scaled cells; enables the direct integration of raster and vector datasets without requiring further rasterization, vectorisation or resampling of data; and, readily facilitates the parallel processing of spatial analyses across high performance and cloud compute environments.

[0031 ] Currently most DGGS implementations define a hierarchical tessellation of the surface of the Earth; however, many of the issues we are increasingly trying to understand, solve and/or manage extend beyond the surface of the Earth to the interior and/or the atmosphere and beyond. To improve our ability to answer these 'Big Data' questions at scale and in near-real time it is necessary to extend DGGS implementations beyond surface implementations and to construct scalable "volumetric" DGGS implementations that leverage and extend the benefits and capabilities of "surface" DGGS.

Advantageous Effects of Invention

[0032] The object of this invention realises an efficient method to construct three- dimensional volumetric model spaces that are scalable from infinitesimally small to global scales. The method of indexing provides a way to map three- dimensional geospatial data to a one-dimensional array representation on a computer data storage device that preserves the local spatial relationships of the data.

[0033] The object of this invention provides an advanced Spatial Data Infrastructure that enables scalable integration and analysis of multiple spatial data sources in three-dimensions using array set theory based on voxel reference identifiers rather than by repeated spatial computations. [0034] The object of this invention extends and enhances the application of OGC compliant "surface" DGGS implementations to enable the resulting Spatial Data Infrastructure to be directly applied to three-dimensional problems, such as city/urban environments, mining and geophysics, oceanography, climate change research and hazard monitoring, to name just a few.

Brief Description of Drawings

Fig.1

[0035] [Fig.1 ] Prior Art Surface Discrete Global Grid System (DGGS) Construction; illustrates the prior art of constructing "surface" surface DGGS reference frames that are conformant with the OGC Discrete Global Grid Systems Abstract

Specification [NPL J] (OGC 15-104r5).

Fig. 2

^00367[Figs. 2A, 2B, and 2C] Equal-Volume Ellipsoidal Shell subdivision of the Earth; illustrates the subdivision of an ellipsoidal model of the Earth (e.g.

WGS84/GRS80) into equal-volume ellipsoidal shells extending from the Earth's centre of mass to an outer radius defined by a positive integer multiple of the Earth's volume.

Fig. 3

[0037] [Figs. 3A and 3B] Example visualisations of a "volumetric" DGGS reference frame constructed using the rHealPIX "surface" DGGS for resolution 0 and resolution 1 respectively; illustrates the construction of a "volumetric" DGGS following the methods of this invention.

Fig. 4

[0038] [Fig. 4] Example resolution 1 subdivision of resolution 0 "volumetric" DGGS voxels; illustrates the volumetric refinement of a set of resolution 0 "volumetric" DGGS voxels following the methods of this invention.

Fig. 5

[0039] [Fig. 5A] Elements of the Spatio-Resolutional Voxel Indexing Schema;

illustrates the key elements required to construct the spatio-resolutional voxel reference identifiers following the methods of this invention (i.e. ellipsoidal Morton Code axes, encoding a spatial location to an ellipsoidal Morton Code of variable lengths from 1 to 64 bits, and the binary representation of the "volumetric" DGGS resolution).

[0040] [Figs. 5B and 5C] Example four-dimensional spatio-resolutional search of a "volumetric" DGGS constructed following the methods of this invention - illustrates the use of 1D space filling curve range searching to return all

"volumetric" DGGS voxels within a three-dimensional search region for a given range of "volumetric" DGGS resolutions.

Fig. 6

[0041 ] [Fig. 6] Typical volumetric Digital Earth Information System. Fig. 7

[0042] [Fig. 7] Typical distributed network of volumetric Digital Earth Information Systems.

Description of Embodiments

[0043] [Fig. 1 ] shows the prior art methodology of constructing a surface Discrete Global Grid System as a hierarchical series of equal area polyhedral tessellations of the surface of the Earth. The method defined by the OGC Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15-104r5) begins with a base polyhedron selected from the list of Platonic Solids (tetrahedron 101 ,

hexahedron/cube 102, octahedron 103, dodecahedron 104 and icosahedron 105). The selected base polyhedron is mapped to the surface model of the Earth (e.g. WGS84 or GRS80) to produce a corresponding initial "surface" DGGS tessellation 106-110. There are three (3) main types of cell shape used to construct a "surface" DGGS (triangles 111 , squares 112 and hexagons 113). The initial "surface" DGGS tessellation is then iteratively refined by a refinement ratio to produce a hierarchical system of discrete global grids; where each cell of the previous level of refinement is divided into an equal number of "child" cells.

Examples of the first two (2) levels of refinement for the corresponding initial "surface" DGGS tessellations (106-110) are shown 114-119 and 120-125

respectively using a refinement ratio of nine (9) (i.e. nine (9) child cells to every parent cell; other valid refinement ratios include four [4], or sixteen [16], child cells to every parent cell). A "surface" DGGS operates through the use of a cell indexing schema (not shown in this figure) where each cell is referenced by a unique identifier; often associated with the individual faces of the base

polyhedron and the specific "surface" DGGS implementation. The construction of a "surface" DGGS that is compliant with the OGC® Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15-104r5) provides a solid base for the construction of a three-dimensional "volumetric" DGGS as described in this invention. 44] [Fig. 2A-2C] illustrates the subdivision of the reference ellipsoid (e.g. WGS84 or GRS80) into equal volume ellipsoidal shells. This is the first step required to extend any "surface" DGGS implementation to the "volumetric" DGGS

implementation of this invention. The reference ellipsoid 2000, and the outer shell of the "volumetric" DGGS 2100 are shown in [Fig. 2A]. The outer shell 2100 represents the same volume above the surface of the reference ellipsoid 2000 as that contained within it. Following the method of Claim 2, [Fig. 2B] and [Fig. 2C] show the subdivision of the ellipsoidal shells shown in [Fig. 2A] using a radial refinement ratio of three (3) (or \j 'surface DGGS refinement ratio = in the examples shown here). In [Fig. 2B] the Reference ellipsoid 2000 is subdivided into three (3) inner shells (2010 and 2020) each representing an equal volume between each shell that is 1 /3 of the volume of the original reference ellipsoid 2000 (here 2000 also represents the surface shell of the subdivided reference ellipsoid) and three (3) outer shells (2100, 2110 and 2120) each representing an equal volume between each shell that is 1 /3 of the volume of the original reference ellipsoid 2000. In [Fig. 2C] each of the shells constructed in [Fig. 2B] (2000, 2010, 2020, 2100, 2110 and 2120) are subdivided into three (3) sub-shells representing 1 /9^th of the volume of the original reference ellipsoid 2000. This results in a new set of ellipsoidal shells 2000, 2001 , 2002, 2010, 2011 , 2012, 2020, 2021 , 2022, 2100, 2101 , 2102, 2110, 2111 , 2112, 2120, 2121 and 2122 respectively. [0045] To demonstrate the construction of a "volumetric" DGGS following the method of this invention, the rHealPIX DGGS will be used as the chosen "surface" DGGS. The rHealPIX DGGS is conformant with the OGC® Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15-104r5). The embodiment of this invention is not limited to the rHealPIX DGGS, any and all flavours/styles of "surface'OGGS which are conformant with the OGC® Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15-104r5) may equally be used to construct "volumetric" DGGS instances following the method of this invention.

[0046] [Fig. 3A] and [Fig. 3B] illustrate how the equal-area cells of a "surface" DGGS (in this case the rHealPIX DGGS) are extended to construct the equal-volume voxels of a "volumetric" DGGS. [Fig. 3A] illustrates the initial tessellation of the "volumetric" DGGS (resolution 0) showing the resolution 0 rHealPIX "surface" DGGS tessellations 301 and 302 on each ellipsoidal shell 2000 and 2100

respectively. The vertices of the corresponding "surface" DGGS cells are joined to form the equal-volume voxels of the "volumetric" DGGS. Each voxel is referenced at its three-dimensional centroid location (e.g. 303). [Fig. 3B] illustrates the resolution 1 refinement of the "volumetric" DGGS showing the resolution 1 rHealPIX "surface" DGGS tessellations 304 - 309 on the

corresponding ellipsoidal shells 2000, 2010, 2020, 2100, 2110 and 2120

respectively. The centroid locations of each voxel 310 are also shown. [Fig. 4] shows the resolution 0 401 and resolution 1 402 refinements of the "volumetric" DGGS for an equatorial resolution 0 rHealPIX DGGS cell (e.g. O, P, Q or R in the rHealPIX indexing schema). The resolution 0 voxel centroid 403 and resolution 1 voxel centroid 404 locations are also shown.

[0047] By translating ellipsoidal coordinates (longitude, latitude, radius) to Earth

Centred Earth Fixed (ECEF) coordinates (x, y, z) it is possible to treat the vertices of "surface" DGGS cells between ellipsoidal shells as points in free space within the same reference frame (i.e. referenced to the Earth's centre of mass rather than to separate ellipsoidal shell reference frames). This simplifies the

specification and computation of voxel vertices, voxel boundaries and voxel centroids. Assigning Identifiers to Voxels

[0048] Just like there are different flavours of "surface" DGGS that are compliant with the OGC Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15- 104r5), there are also many different indexing schemas that may be used to identify and address individual cells. This can introduce added complexity when one proposes to extend specific "surface" DGGS implementations to construct "volumetric" DGGS implementations that also extend the associated "surface" DGGS indexing schemas from a two-dimensional realisation to a three- dimensional one. The indexing schema proposed by this invention avoids this added complexity by defining a simple but effective and completely independent schema that globally references only the position of the voxel centroid and the DGGS resolution level.

[0049] [Fig. 5A] illustrates the use of an ellipsoidal realisation of a Morton space filling curve (Z-curve) 501 to encode the three-dimensional floating point ellipsoidal coordinates of each voxel centroid location into a single integer based spatial reference code (or identifier). By specifying different binary lengths of the Morton code it is possible to encode the voxel centroid location along each Morton code axis to a desired level of precision 502. For example, a Morton code constructed by the methods of this invention with a length of 16 bits 503 will be able to resolve and encode the position of a voxel centroid to a precision of 0.01 1 deg. in longitude, 0.055 deg. in latitude and 245.2 m in radius. Extending the length of the Morton Code to 32 bits 504 will enable it to resolve and encode the position of a voxel centroid to a precision of 1 .67e-7 deg. in longitude, 8.38e-8 deg. in latitude and 0.004 m in radius. For many "volumetric" DGGS

implementations defined by the methods of this invention the preferred embodiment of the voxel indexing schema is for a Morton code which is 64 bits in length 505. This provides a Morton code index with a precision of 3.90e-17 deg. in longitude, 1 .95e-17 deg. in latitude and 8.71 e-13 m in radius (i.e. a three- dimensional spatial precision in the order of picometres to femtometres) and enables the construction of a "volumetric" DGGS that can accurately represent and reference voxels from molecular to global scales; although, for many use cases that level of precision will not be a critical requirement, with precision in the order of millimetres to centimetres (or even courser scales) being adequate.

[0050] To illustrate the construction of an ellipsoidal Morton index by the methods of this invention, for a given location 506, the corresponding Morton index is constructed using Morton code lengths ranging from 1 to 64 bits. For each bit of the Morton code both the Morton code score 507, defined by the method of claim 5, and the location ordinate residuals not yet encoded by the Morton index 508 are shown. With a 64 bit Morton code 509 the spatial residual that is not encoded by the index is 8.673e-19 deg. in longitude, 6.072e-18 deg. in latitude and 3.71 e- 31 m radius.

[0051 ] To link the spatial index of this invention to the "surface" DGGS hierarchy topology, and thus enable easy navigation between parent and child voxels, this invention appends an 8-bit binary representation of the associated "surface" DGGS resolution 510 to the ellipsoidal Morton Index to form a composite four- dimensional spatio-resolutional index that enables the volumetric integration and analysis of data stored at different spatial resolutions without requiring any resampling of the data in question.

[0052] [Fig. 5B] illustrates the use of the composite Morton index to efficiently

perform four-dimensional spatio-resolutional searches of a "volumetric" DGGS constructed by the methods of this invention. A bounding prism/search region 511 is defined. For each vertex of the search region 512 the composite Morton index is encoded 513 for the desired DGGS resolution (in this example

resolutions 0, 1 and 2). This invention leverages the property of space filling curves that enables multi-dimensional spatial queries to be reduced to a one- dimensional query along the space filling curve. For a given search region (e.g. 511 ), all "volumetric" DGGS voxels with a centroid that falls within the search region will have a composite Spatio-Resolutional Morton Index that also falls within the minimum 514 and maximum 515 composite Morton code range defined by the vertices of the search region 516. [Fig. 5B] shows the results of three (3) spatio-resolutional searches of the search region discussed above (i.e. 511 ), for DGGS resolution 0 517, 1 518 and 2 519. [Fig. 5C] shows the results of a spatio- resolutional query for resolution 2 voxels within an extended the search region from the centre of mass of the Earth to the outermost ellipsoidal shell radius 520.

Apparatus and Operating Environment

[0053] The apparatus and operating environment of this invention is illustrated in [Fig. 6] and consists of a volumetric Digital Earth Information System 601 , which can be instantiated on any computational infrastructure with the following hardware devices:

[0054] A central processor unit 602 consisting of at least one (1 ) processor core, but ideally multiple processor cores. The processor unit enables the construction and operation of a "volumetric" DGGS following the methods of this invention;

[0055] A storage device 603 coupled with the processor unit 602 that enables the storage and retrieval of spatial observations and data assigned to the "volumetric" DGGS. The storage device may be a single unit or an array of storage units to enable distributed storage of data within the computational infrastructure that the vDEIS is instantiated on;

[0056] A terminal device, consisting of a monitor and associated graphics hardware

604, which enables visualisation of data stored in the "volumetric" DGGS, with keyboard and pointer devices 605 to enable user interaction with the vDEIS;

[0057] A network adaptor 606 to enable connectivity between the terminal device and one or more processor units distributed across a private network and/or the internet 607. The terminal device may be either local or remotely connected to the storage and compute devices.

[0058] The volumetric Digital Earth Information System 601 described by this

invention may be constructed and operated locally on a single computer, or remotely on a server or cloud/HPC computational infrastructure 607 accessed via a terminal device 604. This makes the application of this invention very flexible and scalable to specific user requirements.

[0059] The ability to distribute data and computation across multiple compute

infrastructures is both a crucial requirement for spatial data infrastructures in the era of "Big Data" and a core capability of DGGS. The volumetric Digital Earth Information System 601 described by this invention leverages this DGGS capability to facilitate distributed and parallel data storage and computation of spatial observations through one, or a number of, "volumetric" DGGS instances constructed using the methods of this invention.

[0060] The voxel indexing schema described by this invention 513 facilitates further flexibility and application of this invention by providing a common schema for the query and transmission of data and information across multiple volumetric Digital Earth Information Systems 601. This is because the indexing schema described by this invention 513 has no direct dependence on the individual indexing schemas of any "surface" DGGS used to construct the "volumetric" DGGS described by this invention. [Fig. 7] illustrates the construction of a distributed network of volumetric Digital Earth Information Systems 701 , each using a different styles of "surface" DGGS 702, across local/private networks and/or the internet that can operate using a single user interface 703 that is based on a common indexing schema 513 and a consistent set of functional algorithms to perform data search, analysis and integration.

Examples

[0061 ] There are no known published examples of the object of this invention in

operation. Pangaea Innovations Pty. Ltd. are currently investigating methods to further develop and bring this technology to market.

Industrial Applicability

[0062] The object of this invention has a very wide industrial applicability as a

scalable spatial data integration engine. Some potential industrial applications include: a. Regional to continental scale geological and geophysical modelling and analyses; b. A near real-time spatial data integrator for Augmented/Artificial Reality technologies; c. Agriculture/Aquaculture environmental monitoring; d. Insurance and hazard risk assessment and mitigation; e. Distributed gaming and model simulation applications in three- dimensional environments; f. Near Real-time monitoring of distribution, logistics, energy network

dynamics in three-dimensions;

Non Patent Literature

[0063] NPL1 : Gibb, R.G., 2016, "The rHealPIX Discrete Global Grid System"

Proceedings of the 9th Symposium of the International Society for Digital Earth (ISDE), Halifax, Nova Scotia, Canada. IOP Conf. Series: Earth and

Environmental Science, 34, 012012. DOI: 10.1088/1755-1315/34/1 /012012

[0064] NPL2: Goodchild, M. F., 1992, "Geographical information science",

International Journal of Geographical Information Systems, 6(1 ): pp 31 -45.

DOI:10.1080/02693799208901893

[0065] NPL3: Goodchild, M. F., 2000, "Discrete global grids for digital earth",

International Conference on Discrete Global Grids. Santa Barbara.

http://www.ncgia.ucsb.edu/globalgrids/papers/goodchild.pdf

[0066] NPL4: Kimerling, A. J., Sahr, K., White, D., & Song, L, 1999, "Comparing

Geometrical Properties of Global Grids", Cartography and Geographic

Information Systems vol. 26, no. 4, pp 271 -88.

DOI:10.1559/152304099782294186

[0067] NPL5: Mahdavi-Amiri, A., Alderson, T., & Samavati, F., 2015, "A Survey of Digital Earth Representation and Visualization", Computers & Graphics, Elsevier Ltd., pp. 95-1 17. uri: http://hdl.handle.net/1880/50407

[0068] NPL6: Peterson, P. R., 2016, "Discrete Global Grid Systems." In The

International Encyclopedia of Geography, edited by Douglas Richardson. Maiden, Oxford: John Wiley and Sons, Ltd.

[0069] NPL7: Purss, M.B.J., Gibb, R., Samavati, F., Peterson, P., Rogers, J.A., Ben, J. and Dow, C, 2017, "Discrete Global Grid Systems Abstract Specification - Topic 21 ", Open Geospatial Consortium, Abstract Specification, 15-104r5. [0070] NPL8: Sahr, K., White, D., Kimerling, A. J., 2003, "Geodesic discrete global grid systems", Cartography and Geographic Information Science, 30 (2), pp 121— 134. DOI:10.1559/15230400310001 1090

[0071 ] NPL9: Yu, J., Wu, L, Li, Z. and Li, X., 2012, "An SDOG-based intrinsic

method for three-dimensional modelling of large-scale spatial objects", Annals of GIS, 18(4), pp. 267-278, DIO: 10.1080/19475683.2012.727865

[0072] NPL10: Zhu, L, Sun, J., Li, C. and Zhang, B., 2014, "SolidEarth: a new

Digital Earth system for the modelling and visualization of the whole Earth space", Front. Earth Sci., 8(4), pp 524-539, DOI: 10.1007/s1 1707-014-0438-7

Claims

A computer-implemented method for extending a two-dimensional hierarchical series of polyhedral tessellations of the Earth's surface (known as "surface" Discrete Global Grid Systems) to define a three- dimensional/volumetric hierarchical series of polyhedral tessellations of the Earth, its atmosphere and outer space extending to an arbitrary radius from the centre of mass of the Earth governed by equal-volume multiples of the volume of the Earth's reference ellipsoid (e.g. WGS84 or GRS80), the method comprising: a. Defining a hierarchical equal-area series of polyhedral tessellations of the Earth's surface reference ellipsoid (e.g. WGS84) in compliance with the Open Geospatial Consortium (OGC) Discrete Global Grid Systems Abstract Specification [NPL:7] (OGC 15-104r5). b. Defining a series of equal volume ellipsoidal shells subdividing the volume of the reference ellipsoid by a fixed number of shells

associated with each successive resolution of the "surface" DGGS. c. Defining the same "surface" DGGS tessellation from a) on each

ellipsoidal shell from b). d. Defining three-dimensional equal-volume voxels by connecting the vertices of each DGGS cell from c) between each adjacent ellipsoidal shell from b). e. Assigning each voxel, a unique hierarchical index comprising an

ellipsoidal Morton Code representation of the voxel centroid location and a binary representation of the "surface" DGGS resolution. The method of [Claim 1 ], wherein the key parameter necessary to define equal volume ellipsoidal shells (i.e. the ellipsoidal semi-major axis a_Sheii ) is determined as follows: a. Define the outer limit for the construction of the "volumetric" DGGS as a positive multiple of the semi-major axis of the Reference Ellipsoid (ad). The "maximum radius multiple" (radius_multiple_max) shall be at least 1 . b. For each successive DGGS resolution the outer semi-major axis for each ellipsoidal shell is determined by the following mathematical formula:

V lref

(Num_shells\

shell J

^■'^■shell

4TT(1 - )

= The flattening ratio of the Reference Ellipsoid

4

Volref ⁼ The Volume of the Reference Ellipsoid i. e. -nalh₀

Num_shells = (pGGS_Refinement_Ratio)^resolution

shell = [1 · ·· (radius _multiple_max * Num_shells)]

rad ius_munip ;_{e max} resolution≥ 0

[Claim 3] The method of [Claim 1 ], wherein the same "surface" DGGS that was

implemented on the reference ellipsoid is implemented on each ellipsoidal shell defined by the method of [Claim 2].

[Claim 4] The method of [Claim 1 ], wherein three dimensional voxels are defined for each ellipsoidal shell, as follows:

a. For each cell of the "surface" DGGS implemented on a given ellipsoidal shell: i. Find the cell vertices, using the functional algorithms

implemented by the "surface" DGGS (these functional algorithms are a requirement of the OGC Discrete Global Grid Systems Abstract Specification [[NPL:7] OGC 15-104r5]); ii. If necessary, convert the cell vertices from angular coordinates (latitude, longitude) to Earth Centred Earth Fixed (ECEF) coordinates using standard geodetic methods; iii. Connect the cell vertices defined in ii. (above) to the ECEF

coordinates of the cell vertices for the corresponding cell defined on the next lower ellipsoidal shell; and, iv. If there are no more ellipsoidal shells below the current

ellipsoidal shell, connect the ECEF coordinates of the cell vertices with the ECEF origin (0,0,0). The method of claim 1 , wherein a unique index is assigned to each voxel, defined by the method of claim 4, as follows: a. Compute the three-dimensional centroid of the voxel in ECEF

coordinates, as follows:

∑(Y oxel_fo_0unci_{ary χ}, Voxel_fo_0unci_aryy, Voxel_fo_0unci_{ary z} )

Voxel_Centroi_d = - n where, n = The number of points along the boundary of the Voxel used

Convert the voxel centroid location from ECEF to ellipsoidal

coordinates using standard geodetic methods.

Encode the ellipsoidal centroid location (longitude, latitude , radius) an integer Morton space filling curve identifier, as follows: i. Construct a Morton code axis subdivision of the longitude

dimension, as follows:

360 360 360 360

Morton _axisi_on— [^■

(2°) ' (2¹) ' (2²) ' ' (2ⁿ) where, n

= The length of the Morton Code in bits (e. g. 32, 64, etc ... ) ii. Construct a Morton code axis subdivision of the latitude

dimension, as follows:

180 180 180 180

Morton_axis_lat = [_ _ _ ... _] where, n

= The length of the Morton Code in bits (e. g. 32, 64, etc ... ) iii. Construct a Morton code axis subdivision of the radial

dimension, as follows:

r₀ r₀ r₀ r₀

Morton_axis_radius = [—,—,—, ... ,—] w/iere, r0 = ^asheii_ma asheii_ma

= equatorial radius of the outermost ellipsoidal shell defined by the method of claim 2: if radius_mui_tipi_emax = 1, a_sheUmax = a₀)

n = 77ie length of the Morton Code in bits (e. g. 32, 64, etc ... ) iv. Adjust the centroid longitude ordinate to be within the range 0.0° to 360.0°, instead of the range -180° to +180° (180.0° of the adjusted longitude range equating to be the origin longitude of the reference ellipsoid [i.e. the Greenwich Meridian 0° E]) v. Adjust the centroid latitude ordinate to be within the range 0.0° to 180.0°, instead of +90° to -90° (0.0° of the adjusted latitude range equating to the North Pole and 180.0° equating to the South Pole) Encode the centroid location as an ellipsoidal Morton Code as follows:

1 . For each bit of the Morton Code (i.e. for each

corresponding element of the Morton axes defined above): a. Initialise the Morton Code bit with the value of 0. b. If the longitude is greater than or equal to the

corresponding bit/element of the Morton_axisi_on array add 4 to the corresponding Morton Code bit. And then subtract the Morton_axisi_on array element value from the longitude ahead of the next iteration - this ensures that the longitude ordinate is represented by only the appropriate bits/elements along the Morton_axisi_on. c. If the latitude is greater than or equal to the corresponding bit/element of the Morton_axisi_at array add 2 to the corresponding Morton Code bit. And then subtract the Morton_axisi_at array element value from the latitude ahead of the next iteration - this ensures that the latitude ordinate is represented by only the appropriate bits/elements along the Morton_axisi_at- d. If the radius is greater than or equal to the corresponding bit/element of the Morton_axiSradius array add 1 to the corresponding Morton Code bit. And then subtract the Morton_axiSradius array element value from the radius ahead of the next iteration - this ensures that the radius ordinate is represented by only the appropriate elements along the Morton_axis_rad,us.

The resulting Morton Code bit will have one of the following values for each bit: e. 0 for no tri-ordinate (i.e. longitude, latitude, or

radius) hits; f. 1 for a hit on radius only; g. 2 for a hit on latitude only; h. 3 for a hit on radius and latitude; i. 4 for a hit on longitude only; j. 5 for a hit on longitude and radius; k. 6 for a hit on longitude and latitude; and,

I. 7 for a hit on longitude, latitude and radius d. Append an 8-bit binary representation of the associated "surface"

DGGS resolution to the tail (right hand) end of the voxel centroid Morton code (e.g. resolution_0 = 00000000, resolution^ = 00000001 , resolution_2 = 00000010, etc . ). A volumetric Digital Earth Information System comprising: a. One, or more, data storage devices; b. At least one processor unit coupled to the data storage device (a), the processor configured to: i. Define a hierarchical series of volumetric polyhedral

tessellations of an ellipsoidal model of the Earth at successively finer resolutions (a "volumetric" DGGS), each voxel representing an equal volume in space; ii. Assign a unique identifier to each voxel following the method of claim 5; iii. Store the unique identifiers for voxels in the data storage device; iv. Store data related to a three-dimensional geospatial location or observation within the Earth, on the Earth's surface, within the Earth's atmosphere, or in outer space extending to an arbitrary radius from the centre of mass of the Earth, defined by an integer multiple of the volume of the Reference ellipsoid model of the Earth, in association with a voxel representing the same geospatial coverage on the storage device; and, v. Search and access data stored within one or a series of voxels on the storage device and present, or transfer, it to the user in a form that can be readily interpreted as is or used in subsequent quantitative and/or qualitative analyses. c. A terminal device configured with: i. A monitor and associated graphics hardware to enable visualisation output from the processor unit; ii. A keyboard (either hardware or software coupled with a touch screen monitor) and a pointer device to enable input to the processor unit; iii. A network adaptor to enable connectivity between the terminal device and one or more processor units distributed across a private network and/or the internet. ] A system of claim 6, wherein multiple "volumetric" DGGS instances, defined by the method of claim 1 , may be constructed as a distributed network of volumetric Digital Earth Information Systems connected via the internet, and/or a private network.