OGC Testbed-16: DGGS and DGGS API Engineering Report

A Discrete Global Grid System (DGGS) represents a spherical partitioning of the Earth’s surface into a grid of cells (or zones) (Wikipedia). The OGC Members approved and maintain an Abstract Specification (AS) that captures the foundational concepts for DGGS (OGC 15-104r5). This Testbed task aims to begin the process to move towards an OGC Implementation Standard for DGGS through the creation of open-source DGGS reference implementations. Testbed-16 represents the initial effort of what is considered a multi-initiatives process.

DGGS offer a new way for geospatial information to be stored, visualized, and analyzed. Based on a partitioning of the Earth’s surface into a spherical grid, DGGS allows geospatial information to be represented in a way that more intuitively reflects relationships between data and the Earth’s surface. With DGGS, providers and consumers of geospatial information can eliminate many of the uncertainties and distortions inherently present with traditional coordinate systems. To fully realize the benefits of DGGS, standard-compliant implementations are required to allow zone-ID management across DGGS with varying structure and alignment.

DGGS presents an opportunity for the geospatial community to implement a representation of Earth that is vastly different from traditional coordinate system based approaches. DGGS has the potential to enable storage, analysis and visualization of geospatial information in a way that more accurately reflects the relationship between data and the Earth. While the OGC DGGS Abstract Specification captures fundamental DGGS concepts, the Testbed-16 DGGS thread initiated work to more concretely demonstrate DGGS in order to drive adoption. This includes advancement towards development of a DGGS reference implementation.

Key questions addressed by the work include:

The participants in the Testbed-16 DGGS thread approached these questions by:

Although DGGS implementations have been deployed in previous testbeds, this was the first time that a DGGS thread has appeared in an OGC Testbed. As a consequence, there were widely differing expectations both across the participants and between the participants and the sponsors. However, once all the participants had gained a common understanding of DGGS, progress was remarkably fast and a lot simpler than was anticipated.

Key highlights from this effort include:

Both server implementations adopted a common strategy of wrapping the chosen DGGS libraries into their systems to present the DGGS Reference System (RS) as a collection of virtual vector layers. Each layer (or item in the collection) corresponds to a level of the DGGS RS grid hierarchy. In both API implementations this allowed the client to select whether it wanted a GeoJSON or a JSON payload. With the GeoJSON payload, the DGGS library created the necessary coordinates for each zone on-the-fly. With the JSON payload, a native DGGS geometry was supplied comprised only of DGGS zone-IDs (the unique identifier associated with each zone ).

For the GeoServer OGC API - Processes implementation this means that all the native GeoServer vector processing functionality, and all the existing GeoServer services - such as OGC WMS and OGC WFS - are exposed to traditional GeoServer clients to use DGGS data. A single GeoServer DGGS wrapper library was created to deliver both H3 and rHEALPix solutions. The no-SQL database ClickHouse was chosen as the GeoServer backend datastore for the DGGS data. This proved very fast for delivering Sentinel 2 data and generating both statistical summaries and index calculations (e.g. NDVI, NDBI & NDWI) on Sentinel 2 data. In the backend ClickHouse datastore, the DGGS zone ID was the primary index for the pixels of Sentinel 2 data.

In the OGC API - Features implementation, a linked data approach was taken using an Resource Description Framework (RDF) datastore. This proved very successful and this Testbed experience is already feeding into the OGC GeoSPARQL Standard roadmap. Vector data for Australia’s Level 1 Statistical Meshblocks and Hydrological Catchments were converted to DGGS RDF data with the DGGS Cell-ID as the predicate (i.e primary index) for the data. The OGC API - Features end-point presented this data as either JSON Linked Data (JSON LD) - using the DGGS zone -ID as the only representation of geometry, or as GeoJSON with the DGGS library generating vector geometries on-the-fly as required.

Participants started with an expectation that the Testbed-16 thread would deliver 'a simple application'. That participants were able to stand up an application delivering native DGGS vector and earth observation capability to users of both traditional Geographic Information Systems (GIS) clients and to native DGGS clients reinforced the versatility of DGGS and underscored its promise to deliver rapid multi-disciplinary analyses. Systems based on DGGS have been shown to have the characteristics of a chameleon that can acquire the character of both vector GIS systems and of earth observation systems. This means that both raster and vector data can be served through either raster or vector or DGGS native services. In its native form, DGGS also offer additional analytical versatility that is not present in either of these.

Testbed participants documented future tasks in the following areas:

All questions regarding this document should be directed to the editor or the contributors:

Contacts

Attention is drawn to the possibility that some of the elements of this document may be the subject of patent rights. The Open Geospatial Consortium shall not be held responsible for identifying any or all such patent rights.

Recipients of this document are requested to submit, with their comments, notification of any relevant patent claims or other intellectual property rights of which they may be aware that might be infringed by any implementation of the standard set forth in this document, and to provide supporting documentation.

The new draft DGGS Abstract Specification (OGC Topic 21 v2.0 / ISO/DIS 19170-1) defines 1.) a DGGS Reference system as 'referencing by zonal identifiers with structured geometry', and 2.) a DGGS as a holistic system comprising:

Most existing DGGS libraries pre-date the new draft Abstract Specification (AS) and as a consequence they only partially fulfil the requirements of the draft AS. Furthermore, they are not structured neatly into modules that follow the above pattern. This has created confusion as to what comprises a DGGS. As such participants in this thread held diverse views of what constituted a DGGS. As the diversity was explored the need for some additional concepts or terms that reflected the different roles that a library might play in DGGS were identified. The participants also identified that some existing libraries supported a single DGGS RS, while others supported whole families of DGGS RSs. An analogy might be the difference between a library that locked in the Universal Transverse Mercator (UTM) Zones 3 projection as compared to a library that supports Transverse Mercator projections. With appropriate configuration parameters the library can support any UTM Zone - including UTM Zone 3. Following the structure outlined above, the participants distinguished between libraries that were DGGS RS providers - and those that were DGGS RS navigators.

Early in the discussion between participants, the question was raised as to whether there was such a thing as DGGS data, and if there were, what would a DGGS datafile look like? For anybody coming into DGGS from a traditional GIS background, and thinking in terms of DGGS as a reference system, the obvious answer was no. For example, concepts such as World Geodetic System 1984 (WGS84) data or UTM 3 data have very little meaning. On the other hand, those participants heavily involved in DGGS thinking were convinced that DGGS data was a very real thing. Thinking about spatial data is typically in terms of its storage, and storage design is intimately bound to the form of the geometry being stored — Shapefiles, HDF5, GeoJSON, LASer (LAS) formatted files — are all solutions to the needs of storing geometries. Traditional reference systems are independent of geometry, and as a consequence the choice of reference systems is independent of the choice of data format.

By contrast, DGGS reference systems define their own geometry. This geometry is implied by the ZoneID. So, within a DGGS a geometry does not need to be explicitly provided as coordinates, and therefore does not need to be encoded in a traditional spatial data format. Instead, DGGS data only needs to record the ZoneId(s) associated with the feature, observation, or raster data zone (cell). This is very much like any other unique identifier associated with a list of attributes, and so DGGS data can be stored in many typical datastores. What makes it spatial data is the presence of the DGGS RS’s navigator functions that perform spatial topology operations using the ZonalID(s) in lieu of coordinate geometry. As a consequence, explicit read/write support for DGGS data in a DGGS library is not required. What is expected is handling of the quantization roles that are defined by the draft OGC DGGS Abstract Specification. These roles tell us the geometric relationship between a record and a ZoneId. They include a single coordinate (the zone centroid), or an area (enclosed by the zone’s boundary), or a tile.

A stock-take of DGGS libraries was performed and the following are the libraries found.

To demonstrate that the DGGS work undertaken was applicable to multiple DGGS, the decision was made to implement two DGGS, and the libraries chosen were H3 and rHEALPix.

H3 only supports one DGGS RS, so there are no further parameter choices to be made. The H3’s reference system is referred to as H3RS.

rHEALPix provides a number of hard-coded references systems, such as for covering spherical vs ellipsoidal earth models. For example, AusPix uses the default hard-coded WGS84 ellipsoidal reference system based on the zero meridian. Since Testbed-16 is developing solutions that may be demonstrated in a number of continents, the decision was made to choose a rotated version of the WGS84 that places all the corners of the initial cube in the sea. This reduces the area of land covered by the most distorted zones. This reference system as TB16Pix.

The code to provide TB16Pix, print the definition and export zone centers and edges for visualization in Google Earth is:

As illustration the corners of the northern and southern zones are shown in the next two figures:

The participants had considerable discussions about the context in which an API that talked DGGS would be used. A couple of issues were central to the discussion:

The alternative contexts discussed were:

There was also discussion about whether to implement an API that conforms to OGC API - Features or OGC API - Process. While OGC API - Features would be useful for allowing people to browse a DGGS and has clarity of meaning, it would also potentially be very chatty in use for analytics (i.e. a potentially very large amount of data, in the form of JSON documents, being transferred either as very large OGC API - Features collection documents, or many transactions of individual OGC API - Features documents for each feature of a collection). By comparison OGC API - Processes is relatively opaque, but fits very well with the DGGS premise that DGGS are particularly suited to Big Data applications where processing should be passed to the data. This is as opposed to the case in which features are passed to the client for processing.

For an API that is a profile of the draft OGC API - Process, the basic resources "are" Processes. Any particular API and API implementation offers particular processes or process types, organizes them in characteristic structures, and supports returning particular representations. The most common representations will be 1) descriptions of the processes, or 2) the result of a process "job" invoked by a request.

What really makes an API, though, is how those resources are organized. The trend with OGC API - Processes (via GET operations) is to group processes according to the primary dataset that they act on, what might be termed the process "affordance". This matches up well with the draft OGC API - Common pattern of dividing resources into collections. Depending on whether a hierarchical organization is supported, each collection is either a collection of collections, or a dataset composed of data elements. The sleight of hand, though, is that a dataset in a Processes API is not a target resource, but a means of organizing processes that are able to act on it.

What makes sense is that a DGGS be thought of as a dataset, with or without data properties, that organizes retrieval, linking, and processing types. The dual nature of a DGGS as both a processing engine and a dataset/collection of geometry objects creates flexibility in the way an OGC API for DGGS could be structured as:

These concepts require more elaboration during subsequent OGC standards development and/or OGC Innovation Program activities to standardize the concept of an OGC API for DGGS.

As a result of these discussions the decision was made to implement DGGS in the context of two distinct APIs.

Acknowledging the conundrum of how to deliver geometry, the following strategy is being used by both API implementations:

These are the operations defined in OGC Topic 21 Part 1 v2.0 and ISO 19170-1 Core for ZoneQuery.

Addressing the Call for Participation (CFP) clause D139 DGGS Enabled Data Services, the OGC API - Features implementations for DGGS datasets should support:

The specifics of these implementations are documented in the next section.

As per the Testbed Call for Participation (CFP), the OGC API - Processes implementation in this task supported "the geographic location to zone-ID(s) and reverse conversion", as well as other eventual additions to allow better comparison and discussion.

The participant discussions and technology implementations conducted during this Testbed activity evaluated the question of "is there a justification, or need, to draft a separate implementation standard to codify OGC APIs for DGGS?". The work done in both D137 and D139 allowed this question to be asked from a number of perspectives.

While the "system" and "data" aspects of a DGGS can be implemented using conventional OGC API - Processes and OGC API - Features schemas, there are some subtle, but very important, differences that separates the operation/access of a DGGS resource from those of a conventional Process or Feature service.

DGGS both simplify the traditional view and add new capability.

The simplification can be summarized as:

New capabilities include:

The Testbed 16 Call for Participation (CFP) alluded to non-specific Use Cases that were to demonstrate two possible 'implementation scenarios' involving the following actors:

Based on the CFP, three initial use cases were posed: Use Case #1 - GPS Location to DGGS Cell…​Use Case #3 - Bushfire Impacts from the "Black Summer" Bushfires in Australia. Participant discussions matured the use cases and at the same time participants gained confidence in what was possible both from an API perspective. This led, from an implementation perspective for the server and client endpoints, a fourth much more ambitious use case Use Case #4 - A DGGS version of a DAPA Use Case making use of one or more Jupyter notebooks. This last use case is the one the participants fnally pursuedin the Testbed.

The fourth use case was developed by analyzing the Testbed-16 DAPA thread Use Cases [5] from a DGGS implementation perspective, identifying a set of common elements, and bringing them together as a single Use Case. This Use Case was intended to demonstrate an alternative DGGS roadmap for the Testbed-16 DAPA thread activities.

The GeoTools library and GeoServer web application had no support for DGGS at the beginning of the testbed. Three DGGS libraries have been evaluated in order to establish a baseline for a generic DGGS Java API.

In particular the following libraries were evaluated:

Table 34 compares how the three libraries represent the basic DGGS concepts and the integration with geographic coordinates and basic geometries is interesting.

Table 35 compares the libraries in terms of functionality provided. The following table provides a set of links to tagged versions of the respective software, with eventual comments when the method in question has unusual behavior.

Finally, given the target of integration with GeoServer, evaluating the suitability of integration with a Java Web Server was important.

As a result of the above tests, and considering the time limits of a Testbed, the OpenEaggr library was dropped from the actual DGGS integration experiments.

The rHEALPix integration development produced functional code, with a couple of significant limitations:

The H3 integration has showed no issues so far, being stable, fast and linearly scalable.

The DGGS subsytem in GeoServer is based on a pluggable extension point called DGGSInstance. The interface exposes the basic operations of a DGGS, allowing upstream code to operate against a DGGS without relying on its implementation details:

The DGGSInstance implementations can be discovered at runtime using the Java Service Provider Interfaces (SPI) mechanism. This allows creating new DGGS integrations, packaging them as JAR files, and having GeoServer discover their presence.

In GeoTools and GeoServer a DataStore is a Java interface " used to access and store geospatial data in a range of vector formats".

A DataStore provides access to a set of feature collections, enabling exploration of their structure (list of attributes and type), as well as to query them, and optionally, to modify their contents. The operation model of the DataStore is heavily influenced by the original WFS 1.0 design, and is still a good match for the existing OGC API - Features.

In GeoServer DataStore instances can be plugged in, discovered and configured at runtime using the SPI mechanism. This allows publishing data over various OGC protocols, such as WFS, WMS, WMTS and the OGC API equivalents.

The first DGGS integration with GeoServer was thus the DGGSGeometryStore, exposing the DGGS zones as features, with an identifier (the ZoneId), a resolution, and a geometry (the boundary).

The store currently takes as a parameter the DGGS to be used. In the future configuration parameters for DGGSs like rHEALPix will be exposed as well. The current rHEALPix implementation exposes the configuration for the TB16Pix RS.

The store then exposes a single feature collection, with the same name as the DGGS, that can be configured as a layer in GeoServer, and then consumed from services such as WMS and WFS.

An example WMS request could be:

https://tb16.geo-solutions.it/geoserver/dggs/wms?STYLES=&LAYERS=dggs%3ATB16-Pix%2Cdggs%3Aworld&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=-106.875,-45.3515625,163.125,89.6484375&WIDTH=768&HEIGHT=384

The following images show maps of both H3 and TB16-Pix from the WMS server:

A WFS request can be used instead:

https://tb16.geo-solutions.it/geoserver/dggs/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=dggs%3ATB16-Pix&maxFeatures=1&outputFormat=application%2Fjson

Which results in:

When using a WMS endpoint, the geometry store dynamically chooses which resolution level to return based on the current scale denominator.

The resolution level can be controlled via the viewparams vendor parameter. This allows passing key values to the data sources. Originally conceived to expand variables in layers sourced from SQL statements, the viewparams mechanism is now available to every store.

In particular, the resOffset view parameter can be used to offset the target resolution. This can be used, for example, to display the DGGS layer multiple times, and viewing the parent/child relationship in a WMS display:

A sample URL is URL-decoded and dissected below:

Static displays of the parent/child geometric relationships in H3 and TB16-Pix follow.

To align with the DAPA tasks, testbed participants decided to expose Sentinel 2 data over the Australian Capital Territory (ACT), resampled over the H3 and TB16-Pix zones, at resolution 11 and higher. For each zone, the following information is available:

To ensure enough data coverage, a bounding box slightly larger than the ACT has been used, in order to ensure enough coverage for the higher level parent zones. Figure 8 shows a reference display of the area.

Even with such a small area, and reduced set of resolutions, the number of records to be stored is significant. In addition to that, to exercise the DGGS and DAPA API at multiple times, snapshots of the area have been taken in 2019 and 2020.

Table 36 & Table 37 report the number of actual zones imported per level (data import is discussed later):

This clarifies the requirements on the data storage, in particular:

A classic relational database can still handle 20 million records at ease. However, with larger areas or deeper resolution sampling, the number would grow very large. The whole world sampled at the highest resolution level of H3 would require storing 597 trillion records.

For the use case, a dedicated On-Line Analytical Processing (OLAP) database is probably a better choice. ClickHouse was the choice to demonstrate this option. The following characteristics make ClickHouse a good match for the use case:

Unlike other recent OLAP database, such as QuestDB, the partitioning key is free to use any column designed, thus supporting splitting the data over both space and time. The following is an example of table creation, using the simplest partitioning engine available in ClickHouse:

Highlights:

This setup allows to quickly resolve queries like the following, extracting all the children of zone P57624 at resolution 10.

For reference, the above structure supported storing all zones in the test area — around 10 million per DGGS type — in roughly 500MB per table. This is a couple of times larger than an equivalent compressed raster storage covering the same area ^[15], at the same resolution. However, there is advantage of a fast SQL aggregation engine built-in, and the freedom to mix columns of different data types.

The following query, computing the maximum value of the B01 over the 8 million zones covering the test area at zoom level 11, ran in 0.015 seconds ^[16]:

As part of the Testbed activities, GeoSolutions imported Sentinel 2 data covering the test area over two time slices, once for H3, and once for TB16-Pix.

The data import was performed using the following steps:

In particular, the Java program for H3 used 16 parallel threads to perform the sampling, loading the 9 million zones in the database in 25 seconds. The generation of all upper levels, by aggregation query in ClickHouse, completed in 1.5 seconds.

The following sample query created all resolution 10 H3 zones, aggregating the values of their direct children at resolution 11, but only when all the children of a given parent zone were available. The query leverages ClickHouse native H3 support functions:

The equivalent rHEALPix query can be built using simple string manipulations (the parent of a zone can be obtained by removing the last character from the zone id) and has a simpler "having" condition, considering a rHEALPix zone always has 9 children:

A second data store was written in GeoServer, reading data from ClickHouse tables satisfying the following structure:

Similar to the geometry DGGS store, this store has no ties to a particular DGGS implementation, and requires the administration to declare which one to use when connecting to a database, in addition to common parameters such as database server host, port, database name, username and password.

The store then makes available to GeoServer all tables matching the above conditions as feature collections, which can be configured as layers and consumed via the classic OGC services.

The figures in Table 38 visually compares H3 and TB16-Pix maps of the NDVI index over the test area, gently zooming in first, and then displaying a lake at the highest resolution. One of the two systems looks more fine-grained than the other. This is actually just an illusion, the issue being that the size of the zones in the two systems are not aligned.

An interesting observation is to see how the store translates the bounding box request from a WMS endpoint into an efficient ClickHouse query. Here is a reference WMS request:

https://tb16.geo-solutions.it/geoserver/wms?STYLES=ndvi&LAYERS=dggs:s2-rpix&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&TRANSPARENT=true&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=148.69720458985,-35.911560058594,149.52117919923,-35.087585449219&WIDTH=600&HEIGHT=600

Table 39 shows the resulting rendering, and on the side, a map with the zone identifiers:

The resulting database query looks as follows:

In other words, the system recognizes some of the zones are fully covered by their parent, and as a result, a query against the parent prefix identifier is done, selecting only the zones at the desired resolution.

Using the viewparams=resoffset:1 vendor parameter to retrieve zones at resolution level 7 in the same area does not change the structure of the generated query, the system still recognizes P67307, P67307, P67304 and P67303 are covering most of the map. Therefore the "like" portion of the query is unchanged. However, the "in" portion of it has a different (and longer) list of zones partially intersecting the bounding box.

This approach helps keeping the queries small, which is particularly important when trying to render a large number of zones using the resolution offset vendor parameter.

Once the Sentinel 2 datasets have been loaded in GeoServer, it’s possible to display false color maps of a given area. The following maps show a constrast stretched composition of B04,B03,B02 bands over the whole ACT area:

Focusing on Canberra, the following three maps show different views of the same area:

The ogcapi-dggs module implements the process-oriented DGGS API described in the DGGS process API chapter.

The API is implemented using the GeoServer framework for OGC APIs, shared with OGC API - Features, and Styles and Tiles.

The implementation provides a HTML and JSON representation of most resources, with a few additional formats when required, or there is an opportunity for reuse of existing encoding support.

The following pictures provide a walk-through of the API resources available at the GeoSolutions Testbed-16 server: https://tb16.geo-solutions.it/geoserver/dggs/ogc/dggs

In particular, the pictures show the basics of the API as well as an exploration of the TB16-Pix collection, thereby providing a purely geometrical description of the TB16-Pix DGGS (backed by the ÐGGSGeometryStore implementation).

What follows is a screenshot of the s2-rpix zone listing, containing Sentinel 2 data over the ACT, in September 2020:

The s2-rpix and s2-h3 collections are interesting, as they are time enabled. Time management is configured in the UI as follows:

This affects both WMS instances and the DGGS API. In WMS the capabilities document reports the available times, and defaults to the September 2020 time slice.

At the DGGS API level, the time is exposed in the collection description:

All the API requests in addition can receive a datetime parameter specifying the desired time slice. If not specified, the API assumes the configured default.

The following request extracts a single zone, without a datetime specified, resulting in the September 2020 instance to be extracted:

https://tb16.geo-solutions.it/geoserver/ogc/dggs/collections/s2-h3/zone/?zone_id=86be0d207ffffff&f=application/dggs%2Bjson

Providing a specific datetime allows to retrieve a different time slice:

All API endpoints returning zones accept a datetime query parameter, allowing full usage of temporal datasets.

The following example retrieves resolution 6 zones matching the CRS84 point 149,-35 during September 2019:

https://tb16.geo-solutions.it/geoserver/ogc/dggs/collections/s2-h3/point?point=149,-35&resolution=6&datetime=2019-09&f=application/dggs%2Bjson

Similarly, a polygon mapping operation can be performed, mapping the CRS84 ellipsoidal triangle to a set of resolution 5 zones, for September 2019:

https://tb16.geo-solutions.it/geoserver/ogc/dggs/collections/s2-h3/polygon?polygon=POLYGON149%20%20-36,%20149%20-35,%20150%20-35,%20149%20-36,%20&resolution=5&datetime=2019-09&f=application/dggs%2Bjson