The rHEALPix Discrete Global Grid System: considerations for Canada

The geographic grid is a latitude–longitude parametrization of the sphere or ellipsoid (Amiri et al. 2015) that was originally designed for repeatable navigation and naturally suited to vector data, because the coordinate system references a continuous point field (Sahr et al. 2003; OGC 2017b). To create the quadrilateral cells that make up the geographic grid, equal angle steps are taken along the latitude and longitude axes.

Although the geographic grid has many advantages, it contains inadequacies that make it impractical for many applications. For example, cell centroids are not equal distances apart, cells do not have uniform adjacency, and cells surrounding the poles are actually triangles not quadrilaterals. These issues have been highlighted in several applications, including climate modelling (Randall et al. 2002), numerical weather prediction (Carfora 2007), and oceanography (Hoggarth and Taylor n.d.). Furthermore, grid cells do not have equal area (they gradually decrease in size from the equator to the pole), which complicates analyses such as sampling because cells do not have equal probabilities of being selected (White et al. 1992; OGC 2017b). Finally and perhaps most importantly, with the advent of big data, global datasets are now heterogeneous and high in volume. As a result, combining datasets from different sources (such as those having different sensors or coordinate reference systems) on the geographic grid is computationally expensive and time consuming (OGC 2016).

On the other hand, a DGGS has many advantages that helps solve some of the problems faced by the geographic grid. Firstly, with regards to big data, a DGGS can easily manage, store, and rapidly assemble large volumes of heterogeneous geospatial data, which enables efficient exploration, mining, and visualisation. Secondly, a DGGS can easily combine vector and raster data in the same framework, thereby overcoming some of the difficulties faced by GIS approaches such as the “vector–raster divide” and use of map projections (OGC 2017c). Finally, by adopting a grid of equal area cells, spatial analysis can be replicated regardless of location or resolution, and information collected at a given location can be (i) referenced to the explicit area of its associated cell, (ii) integrated with other cell values, and (iii) provide statistically valid summaries from a sample of cells (OGC 2017b).

There have been many DGGSs devised since the 1980s, and new interesting examples are still being created such as the PlanetRisk DGGS (Sahr et al. 2015). A good summary of DGGSs up to 2003 can be found in (Sahr et al. 2003). Although there are a great variety of DGGSs, it should be noted that a single DGGS cannot currently provide the solution to all modern day geospatial applications. Having said that, they have been used separately for line simplification, optimal path determination, statistically valid sampling, and indexing geospatial databases (Sahr et al. 2003). Furthermore, they can be optimized to meet many global geospatial needs including processing, analysis, visualization, and modelling (Stefanakis 2016; OGC 2017a).

In January 2016, the OGC — an international organization that provides standards to the global community — issued a press release seeking comments on a DGGS standard that may replace those based on geographic coordinate systems (OGC 2016). In October 2017, the OGC officially adopted the DGGS as a new standard (OGC 2017c). This was an important announcement that could revolutionise the way geospatial information is used in the future. In addition, the OGC has set up a DGGS Standards Working Group (OGC 2017a) with the aim of increasing awareness of DGGS advantages and to explore a standard that enables interoperability between DGGSs. This has already prompted research on methods to aid index conversion between different DGGSs (Amiri et al. 2015) with more research likely to follow.

The new OGC DGGS standard entitled “Topic 21: Discrete Global Grid Systems Abstract Specification” includes important information on definitions, conventions, functional algorithms, the core data model, and more (OGC 2017b). One section in particular, namely Global Grid Taxonomy, highlights the possible ways to create a global grid, only a few of which are classified as suitable under the abstract specification. Currently, three prominent options are the Snyder projection, rHEALPix projection, and Quaternary Triangular Mesh (QTM). Notably, each is based on a different cell shape. Figure 1 shows examples of DGGSs that use these methods. Interestingly, PYXIS Innovation (PYXIS Innovation 2017) — a Canadian company that has gained widespread attention for creating a new DGGS platform to analyse geospatial data — uses the ISEA3H (Icosahedral Snyder Equal Area Aperture 3 Hexagonal) DGGS (Sahr et al. 2003), whereas people at Landcare Research New Zealand created the rHEALPix DGGS (Gibb et al. 2016). This reiterates the fact that no single DGGS is optimal for all situations or preferred by all users.

With regards to current software on the market for geospatial applications, PYXIS is leading the way. Its Digital Earth web based application makes it is easy to explore, combine, and analyse different datasets by adopting DGGS technology (PYXIS Innovation 2017). However, like any DGGS, the ISEA3H DGGS consists of specific design choices that mean it is not suited to every application. For example, it uses a hexagonal grid which is less familiar to users and less adaptable to Cartesian coordinate systems and hardware devices than traditional square grids (Amiri et al. 2015). In comparison to the hexagonal grid, quadrilateral grids mirror the geographic grid more closely. An emerging DGGS that adopts the quadrilateral grid approach is the rHEALPix DGGS.

Before delving more into the rHEALPix DGGS, a few key terms should be clarified, with the first being base polyhedron. This refers to one of the five platonic solids that can tessellate a sphere or ellipsoid into equal area cells (Fig. 2).

The second term is refinement (sometimes called aperture). This refers to the process of subdividing a cell into a certain number of smaller cells. For example, if a planar square cell is subdivided into four square cells of equal area, the refinement is called one-to-four, nine square cells is called one-to-nine, and so on (Fig. 3).

The third term is resolution. This refers to the number of times the cell subdivision process has occurred. For example, assuming one-to-nine refinement, resolution one for a square has nine cells and resolution two has 81 cells (Fig. 4).

With regards to cells, key terms include parent cell, child cell, aligned, congruent, and adjacency. A parent cell is a cell at resolution, k, and a child cell is a cell at resolution, k+1. Aligned means the centroid of the parent cell is also the centroid of a child cell. For example, this is true of a 3×3 (one-to-nine refinement) square cell but not for a 2×2 (one-to-four refinement) square cell (Fig. 5).

Congruent means the parent cell can be subdivided in such a way that child cells are enclosed by the parent cell. For example, squares and triangles are congruent but hexagons are not (Fig. 6).

Adjacency refers to how a cell is connected to its neighbours. For example, hexagons exhibit uniform adjacency because they have six equal edge neighbours. Conversely, squares and triangles exhibit non-uniform adjacency because they have both edge and vertex neighbours (Sahr et al. 2003; Fig. 7).

In order for the rHEALPix DGGS to be more aligned with geospatial applications, the team made a few important amendments, including

Because reference ellipsoids (such as GRS80 and WGS84) are used in geodesy rather than spheres, the rHEALPix DGGS is more aligned with geospatial requirements and terrestrial use than its predecessor (Gibb 2016). Furthermore, by adopting a platonic solid as the base polyhedron, it is aligned with the criteria and design choices outlined in previous work on DGGS requirements for a digital earth framework (Sahr et al. 2003). Finally, by altering the planar projection to an alignment of square grids, it is easier to understand and link to similar grids used in remote sensing and environmental modelling (Gibb 2016). For detailed information on the rHEALPix DGGS, including the definition, cell IDs, cell shapes, projection equations, distortions, and mathematical underpinnings (see also Calabretta and Roukema 2007), readers are referred to Gibb et al. 2016.

The fundamental features of the rHEALPix DGGS that set it apart from other DGGSs include (Gibb 2016) the following:

Two features that distinguish the rHEALPix DGGS from the popular the ISEA3H DGGS are congruent cell structure and isolatitude distribution of cells. Firstly, congruency is an important property because it makes hierarchy-based cell indexing a trivial process, which enables efficient handling of large volumes of data (Amiri et al. 2013; Amiri et al. 2015). Secondly, isolatititude distribution of cells is a rare quality in DGGSs. As this feature is essential for fast computation of spherical harmonics, the rHEALPix DGGS is a good choice for applications that require harmonic analysis (Gibb et al. 2016). Therefore, geospatial applications involving tidal analysis or image processing, for example, may find the rHEALPix DGGS well suited. Overall, the rHEALPix DGGS contains advantages both from its fundamental design and as a by-product of being a quadrilateral cell based grid. Furthermore, it contains features that other DGGSs simply do not have. As a result, the rHEALPix DGGS may be the ideal choice for new geospatial applications in the future.

It should be mentioned that the rHEALPix DGGS is not perfect for all applications and has limitations of use just like any other DGGS. The main disadvantages are lack of uniform adjacency and varying cell shape. The former, which is common in both quadrilateral and triangle cell based DGGSs, means distances to neighbouring cell centroids are not equal. This causes problems for dynamical systems where functions are related to intercell distances, as well as in applications involving discrete simulations (Sahr et al. 2003; Gregory et al. 2008). The latter is a result of combining four triangular Collignon projections into a single square (Fig. 8). Although all cells are quadrilaterals in the plane, some change shape when projected to the ellipsoid and four different cell shapes occur: the quad cell (quadrilateral), the dart cell (triangular), skew-quad cell (quadrilateral), and cap cell (circular) (Gibb et al. 2016). This is reminiscent of the geographic grid (triangular cells at the pole) and may lead to issues involving conversion algorithms between DGGSs or representation of geometries. Interestingly, however, the rHEALPix DGGS can be rotated, which means some cell shapes can be avoided in a given target area.

In 2016, the rHEALPix DGGS was updated in response to the proposed OGC standard. The work included two main features. Firstly, it showed that cell geometries and cell IDs agree with the OGC standard. Secondly, it defined a process to determine cell adjacency and DE-9IM topological relationships using unique cell IDs rather than geodetic coordinates (Gibb 2016). This update was important work that promoted the usability and benefits of the rHEALPix DGGS, as well as its ability to be compliant with the proposed OGC standard.

In addition, the rHEALPix DGGS has been implemented in Python, and both the rHEALPix and HEALPix projections have been implemented in the Proj.4 Cartographic Library (Gibb et al. 2016). The Proj.4 library allows users to convert geodetic coordinates into and out of both projections using a command line window such as the OSGeo shell. Furthermore, it establishes the ground work for meeting the OGC requirement of having operations that allow data to be assigned to and retrieved from cells (Gibb 2016).

Research yielded one more example of implementation. A document prepared by the Spatial Data on the Web Working Group (a joint World Wide Web Consortium (W3C) – OGC project) described how earth observation data can be represented using the W3C Resource Description Framework (RDF) Data Cube and the rHEALPix DGGS. In particular, there is proof of concept that shows how satellite imagery can be retrieved from the Data Cube using a SPARQL query system (W3C 2017).

Currently, no further work or application utilising the rHEALPix DGGS can be found. With regards to modelling, visualisation, and analysis of geographic entities, PYXIS have shown that DGGSs are a powerful tool. Could the rHEALPix DGGS provide new insights or applications in this field? As of yet, no research has been carried out. Future work should aim to determine how well the rHEALPix DGGS is suited to handling geographic entities such as points, lines, and polygons. Furthermore, due to the lack of open source algorithms that model geospatial entities on DGGSs, algorithms should be made available to the public to promote DGGS usability and reproducible science.

The polar region consists of the north polar region (ϕ > ∼41.9°) and south polar region (ϕ < ∼−41.9°). Because they are isomorphic (Gibb et al. 2016) and Canada (almost entirely) lies in the north polar region, the discussion will be based there almost entirely, because Canada’s most southerly location is in Ontario just below the equatorial and polar region border (e.g., Pelee Island falls in the equatorial region). This is an important consideration if the area of interest is southern Ontario. Essentially though, all of Canada falls within the north polar region, so equatorial region characteristics are not applicable.

The north polar region is created by rotating four triangular Collignon projections from the HEALPix projection into a single square. As a result, parallels project to straight lines (and appear as concentric squares), and meridians project to straight lines heading poleward (Fig. 8b). This representation is quite difficult to understand, and visually determining projection coordinates from geodetic coordinates is not an easy task. Although the mathematical formulas have been implemented in the Proj.4 library to carry out the conversion, this tool should be used with caution. Results can sometimes be unexpected and unfathomable without knowledge of the projection. Conversely, a square grid in the plane is projected to an ellipsoidal grid containing a variety of cell shapes. For refinement one-to-nine, there are three different cell shapes: the cap cell (circular), dart cell (triangular), and skew quad cell (quadrilateral) (Fig. 10). Note, if refinement is even, then there are no cap cells (Gibb et al. 2016).

The cap cell covers the pole and is the only circular cell. Its edge is a line of constant latitude with the cell centroid located at the pole. This cell is not a major concern for Canada, because as resolution increases, the parallel defining the cell boundary converges towards the pole (as the cell area decreases and the centroid remains at the pole). Although Canada’s northern boundary reaches high into the arctic, the likelihood of encountering the cap cell decreases with every increase in resolution, as shown in Fig. 11.

Triangular dart cells are formed where the edges of each of the four Collignon projections meet. The southern edge of a dart cell is a line of constant latitude and north–east and north–west edges converge poleward to a point (Gibb et al. 2016). Furthermore, dart cells are aligned along four meridians and bounded by two parallels: the parallel that defines the cap cell and the parallel that defines the equatorial and polar region border. This divides the north polar region into four quadrants spanning 90°. Because Canada has a large east–west extent that approaches 90° in longitude, dart cells are unavoidable when considering the whole country (although they can be positioned at the east and west edges of Canada, see Fig. 12). However, if smaller features are of interest such as the province of New Brunswick, dart cells can be avoided entirely.

The remaining cells in the north polar region are skew-quad cells. Skew-quad cells have four edges with the north and south edges aligned with parallels. In addition and only for skew quad cells, cell orientation changes with longitude through a quadrant (Fig. 12). Consequently, cells initially having a north–west slant gradually rotate clockwise to a north–east slant. Similar to dart cells, this effect is unavoidable if considering the whole country but can be reduced over smaller areas.

It has been established that a country the size of Canada in the polar region will have issues with varying cell shape and cell orientation. However, smaller areas can escape some of these issues, but how is this achieved? The answer lies in the versatility of the rHEALPix projection.

An important advantage of the rHEALPix projection is its versatility, which means it can be rotated and altered to meet user needs. To tailor the rHEALPix projection to Canada, or any country in fact, it is possible to alter the rHEALPix projection in two ways. Firstly, the north polar square and south polar square can each be positioned in four different locations. Secondly, the projection image of the ellipsoid can be rotated by shifting the position of the prime meridian (0° longitude), see Gibb 2016.

By rearranging the north polar square, n, and south polar square, s, into four different locations, it is possible to create 16 different views of the projection. For example, Fig. 13 shows the (n,s) = (0,0) and (n,s) = (1,3) positions. Choosing a view is not a difficult task because the mathematics behind the projection enables the positions to be easily defined. The main advantage of this is that it allows the projection to be rearranged to enable a better view of a chosen area (Gibb 2016).

The process of shifting the prime meridian simply involves adding the desired offset to all longitudinal values. This results in rotating the grid cells longitudinally around the ellipsoid. With regards to the polar region, this has two key benefits. Firstly, dart cells can be avoided because the four meridians on which they lie can be selected arbitrarily. Although the four meridians will always have a 90° separation, the positions of the meridians can be rotated to any location. Consequently, dart cells can be avoided for any area with a longitudinal extent less than approximately 90°. Because Canada has a longitudinal extent of approximately 90°, dart cells can almost be avoided by shifting the prime meridian to 50° W. This aligns the dart cells with the 50° W and 140° W meridians, as shown in Fig. 14.

Secondly, skew quad cells with angled orientations can be avoided or even exploited. This is an interesting feature because it provides some flexibility to choose cell orientation over a given area. For example, if an area is north–south aligned such as New Brunswick, skew quad cells having the same orientation can be used and a grid resembling the quad cell grid of the equatorial region can be achieved. Figure 15 shows a skew quad cell grid over New Brunswick at various resolutions where the prime meridian has been shifted to 21° W. However, if an area has an angled orientation such as Cape Breton, Nova Scotia, the grid can be rotated so that angled skew quad cells with similar orientation are used instead, see Fig. 16 where the prime meridian has been shifted to 50° W. Notably, this characteristic only applies for a specific range of orientations. The most angled skew quad cells are adjacent to dart cells and have orientations of approximately N 30° E and N 30° W. Therefore, using angled skew quad cells is most suited to areas having an orientation in this 60° range.

With regards to the equatorial region, shifting the prime meridian rotates the grid so that quad cells can be more aligned with a given area. However, because quad cells are always aligned with parallels and meridians, shifting the prime meridian does not affect cell orientation.