Encoding and operation scheme for the rhombic triacontahedron aperture 4 hexagonal discrete global grid system

ABSTRACT

Discrete Global Grid System (DGGS) is a hierarchical structure with seamless global coverage, supporting processing and analysis of heterogeneous geospatial data. Compared with the most commonly used icosahedron, the rhombic triacontahedron (RT) approximates Earth better and can improve the accuracy of data modeling and expression. However, current RT DGGS studies have adopted a single-resolution integer coordinate scheme within a single surface that cannot achieve multilevel grid analysis. To this end, three adjacent surfaces are combined and a three-axis integer coordinate system is established to reduce the number of times crossing the surface. Then, the principle for partitioning aperture 4 hexagonal grids and a quadtree hierarchical encoding scheme are proposed. Further, this study establishes a fast transformation between the code and three-axis coordinates, implements hierarchical and neighbor code operations, and designs a transformation between the code and geographic coordinates. Experimental results indicate that, compared with the icosahedral hexagonal DGGS, the proposed scheme has a significant efficiency advantage in the transformation between the code and geographic coordinates, while the neighbor code operation efficiency can reach 40.97 times that of HLQT and 4.35 times that of HHUT. The proposed scheme provides a more suitable framework for organizing and analyzing global geospatial data.

KEYWORDS:

• Discrete global grid systems; rhombic triacontahedron; hexagon; encoding and operation

1. Introduction

Geospatial data are currently characterized by increasing widespread spatial coverage, large data volumes, various types and formats, and complex data structures. The data processing capability of traditional geographic information systems (GIS) cannot meet this demand (Goodchild 2018; Lewis et al. 2017; Mahdavi-Amiri, Alderson, and Samavati 2015). Therefore, there is an urgent need to construct new models to support global and massive geospatial data fusion and processing. The Discrete Global Grid System (DGGS) partitions Earth into a seamless multi-scale hierarchy, similar to an electronic table covering Earth (Guo, Goodchild, and Annoni 2020; Peterson 2006; Sahr, White, and Jon Kimerling 2003). This is a highly promising framework for organizing and managing spatiotemporal data and has recently been widely applied in wildfire modeling (Robertson et al. 2020), terrain analysis (Li, McGrath, and Stefanakis 2021), flood risk modeling (Chaudhuri, Gray, and Robertson 2021; Fichtner et al. 2023), hate crime studies (Jendryke and McClure 2021), intelligent geospatial maritime risk analytics (Rawson, Sabeur, and Brito 2021), ocean tidal information extraction (Wang et al. 2023a), and ocean surface current estimation (Wang et al. 2023b).

Traditional latitude and longitudinal grids are highly distorted in polar regions, with a large amount of redundancy when used for data organization. The polyhedral DGGS uses specific methods to recursively partition the surface of a polyhedron, followed by projecting it onto Earth’s surface to form globally continuous and uniformly regular hierarchical grids (Mahdavi-Amiri, Alderson, and Samavati 2015; Sahr, White, and Jon Kimerling 2003), rendering it a popular research and application topic. Five platonic solids are commonly used as the base polyhedron to construct DGGS because each surface is a regular polygon of the same size, which can be simply partitioned into regular cells. In general, the greater the number of surfaces in the platonic solid, the less the distortion introduced when projecting between the polyhedral surface and the corresponding sphere (Sahr, White, and Jon Kimerling 2003; White et al. 1998). Therefore, the icosahedron is the most common choice for base platonic solids.

The cell types for grid partitioning are generally divided into triangles, quadrangles, and hexagons, among which hexagons have more ideal geometric properties, such as adjacency consistency, the highest spatial sampling rate and angular resolution (Bousquin 2021; Mahdavi-Amiri, Harrison, and Samavati 2015; Mechenich and Žliobaitė 2023; Yao et al. 2023). The aperture refers to the ratio of the areas of cells at the adjacent level. There are three types of apertures for hexagonal grids: 3, 4, and 7. Aperture 3 grids can provide smoothest transition between adjacent levels; however, the direction of cells at the adjacent level experiences a 30° periodic rotation (Ben et al. 2018; Mahdavi Amiri, Alderson, and Samavati 2019; Peterson 2006; Sahr 2008). Aperture 7 grids have the advantage of fitting the child cells most closely to the area covered by the parent cell; however, the area of cells at the adjacent level has a huge difference while the direction also experiences a rotation of 19.1° (Bousquin 2021; Uber Technologies 2023). The direction of aperture 4 grids does not change at the adjacent level, which is convenient for establishing a multi-resolution structure. Tong et al. (2013) combined two aperture 4 partitioning methods to design the hexagon quadrature balanced structure (HQBS). Wang et al. (2021) designed the hexagon lattice quad tree (HLQT) using complex numbers, characterized by a higher operational efficiency than that of HQBS. Zhao et al. (2022) proposed the optimized hexagonal quadtree structure (OHQS), achieving code operation based on the $\mathrm{ijk}$ coordinate system. Zhou et al. (2023) proposed the hexagon hierarchy on uniform tiles (HHUT), which simplified the processing method when crossing tiles. All the above studies use the icosahedral hierarchical encoding scheme; moreover, common open-source libraries, such as H3, DGGRID, OpenEAGGR, and Geogrid, select the icosahedron to construct DGGS (Li and Stefanakis 2020).

Nevertheless, the icosahedron has significant limitations (Liang et al. 2022). The icosahedral surface is still very different from the surface of Earth, resulting in great distortion and low precision of grid cells. In addition, its triangular surface is incompatible with the commonly used row–column indexing and rectangle data organization, thus complicating the arithmetic and reducing the data processing efficiency. To address these problems, recent studies have attempted to use Catalan solids with more surfaces to construct grid systems with better geometric properties, thereby revealing new design choices for DGGS. Hall et al. (2020) used the disdyakis triacontahedron (DT) to construct a DGGS, but the surfaces of the DT are all irregular triangles, which is only suitable for partitioning of triangular grids while having complex hierarchical relationships. Moreover, the frequent cross-surface operations between 120 triangular surfaces greatly decrease the efficiency of neighborhood queries for grid cells. Bernardin et al. (2011) developed Crusta as an alternative virtual globe based on the rhombic triacontahedron (RT), allowing users to easily visualize custom imagery and topography. Liang et al. (2022) constructed RT hexagonal grids and established a natural row–column coordinate system to index hexagonal cells within each rhombic surface. Wang et al. (2023) utilized the RT slice and dice equal-area aperture 4 hexagonal grid (RTSDEA4H) to achieve rapid grid generation of irregular areas. Considering the complexity of cross-surface operations and the shape of polyhedral surfaces, RT becomes the ideal choice.

RT is composed of 30 congruent ‘golden’ rhombic surfaces (Weisstein 2023), with 32 vertices and 60 edges (Figure 1(a)). The ratio of the long diagonal to the short diagonal on each surface is the golden ratio. Compared to the icosahedron (Figure 1(b)), RT has the following geometric advantages:

1. More number of surfaces results in a higher degree of fitting with Earth. Therefore, the grid cells achieve more optimized geometric characteristics, facilitating global geospatial data integration and mode calculation of the Earth system. For example, compared to icosahedral grid cells, the standard deviation of compactness decreases by 39.6%; the mean and standard deviation of the length deviation decrease by 44% and 32.8%, respectively; the mean and standard deviation of the angle deviation decrease by 29.1% and 25.6%, respectively (Liang et al. 2022).

2. Each surface is a rhombus of the same size and suitable for partitioning using various cells. A set of adjacent edges of a rhombic surface is a natural two-dimensional coordinate system that can easily describe the location of cells, compatible with the rectangle data organization, thus beneficial for data storage and parallel computing. Furthermore, four-side surfaces lead themselves to elegant and efficient quadtree approaches that are easily adapted to graphics hardware and easier to understand for non-professional users as well (Bernardin et al. 2011), which is beneficial for further applications of DGGS.

Figure 1. Comparison of the rhombic triacontahedron and icosahedron. (a) the RT, and (b) the icosahedron. The figures are visualized in Wolfram Mathematica 13.2.1.

However, all the existing hexagonal encoding schemes for RT grid systems adopt a single-resolution integer coordinate scheme that cannot support continuous multi-scale grid analysis. Furthermore, their studies are limited to one single rhombic surface. To this end, this study aims to design an efficient aperture 4 hexagonal hierarchical encoding scheme based on the RT. The organization of this manuscript is as follows. Section 2 briefly describes the basic idea of the proposed scheme. Section 3 describes the methods used in this study, including constructing the combined structure and description of cells’ location, quadtree hierarchical partitioning and encoding, the transformation between the code and three-axis integer coordinates as well as that between the code and geographic coordinates, and implementation of code operations. Section 4 compares the efficiency of the proposed scheme with that of existing approaches and discusses the reasons for the outperformance of the proposed scheme. Section 5 summarizes the key findings and their broader implications.

2. Basic idea

This study chooses the RT as the base polyhedron to construct a new DGGS. Therefore, all components of the DGGS must be revisited and potentially modified (Hall et al. 2020). The basic idea is illustrated in Figure 2.

Figure 2. Schematic diagram of the research methodology.

First, to simplify the frequent cross-surface operations, three adjacent rhombic surfaces are combined into a structure and the three-axis integer coordinate system is established within each combined structure. Second, a new hierarchical hexagon partitioning and quadtree encoding scheme is designed based on the RT combined structure, which is considered to be the simple and efficient method for the subdivision of rhombic surfaces (Bernardin et al. 2011). Third, a mutual transformation between the code and three-axis integer coordinates is constructed to integrate the advantages of the hierarchical code and integer coordinate code. The hierarchical operations are then achieved directly through the code; the neighbor operations are efficiently implemented by the three-axis integer coordinate when the cells are located in the same combined structure and by replacing Morton digits of the code of cells when crossing adjacent combined structures, which avoids rotation and translation in traditional cross-surface operations. Last, using polyhedral projection and the above methods, a precise transformation between the geographic coordinates and code is achieved.

3. Methods

3.1. RT combined structure and cell location description

For hexagonal grids, there are cross-surface cells at the border of the RT rhombic surface, unlike triangular and quadrangular grids; simultaneously, the space coverage of a single rhombic surface is small, frequently resulting in a span of multiple rhombic surfaces when processing geospatial data for a large area. The complexity and low efficiency of cross-surface computing have become the key factors limiting the application of RT DGGS. Therefore, three adjacent rhombic surfaces are combined into a structure, as shown in Figure 3. RT consists of 10 identical combined structures numbered 0–9 in the sequence.

Figure 3. Combined structure of the RT.

The common vertex of the three rhombic surfaces in the combined structure is considered as the origin, followed by establishing a three-axis integer coordinate system, IJK, along the edges of the rhombic surfaces. The negative coordinate components are corrected to 0 for simpler processing. As shown in Figure 4, the spatial position of all cells within a combined structure can be uniformly described without considering the cross-surface problem between the three rhombic surfaces. Thus, the frequent cross-surface problems of 30 rhombic surfaces in the RT can be simplified into operations between 10 combined structures, and this operation is described in Section 3.4.2; the cells within the same combined structure can perform neighbor operations through simple addition and subtraction of integer coordinates, which improves the processing ability of RT DGGS for large-area and large-volume geospatial data.

Figure 4. Diagram of the three-axis integer coordinate system.

3.2. Hierarchical partitioning and encoding

The three-axis integer coordinate has only a single resolution, and establishing a hierarchical structure that reflects the affiliation of cells is impossible. Therefore, this study designs a hierarchical partitioning and encoding based on the three-axis integer coordinate system.

To cover Earth in a regular manner, this study stipulates that the cells overspread the right boundary of the combined structure at each level. In other words, the cells on the right four boundaries of the combined structure belong to the same structure, whereas the cells on the left four boundaries belong to the adjacent structure. Figure 5(a) shows the initial cells in the combined structure, where each initial cell corresponds to a rhombic surface. The white hexagon corresponds to the upper rhombic surface, numbered 0; the gray pentagon corresponds to the middle rhombic surface, numbered 1; and the blue hexagon corresponds to the lower rhombic surface, numbered 2. Different aperture 4 hierarchical partitioning and encoding methods are adopted according to the number of initial cells. As shown in Figure 5(b), the children at different positions are represented by 0–3. Children of white and blue initial cells adopt hexagonal quadtree partitioning (Zhang, Nie, and Chen 2017) and the Morton encoding scheme, ensuring three adjacent boundary children. In contrast, children of gray initial cells adopt another hierarchical partitioning and encoding scheme, such that the child numbered 3 is exactly located outside the current cell. Additionally, the two initial pentagon cells at the south and north vertices of RT are numbered 3, which belong to the No. 2 and No. 7 combined structures, respectively. They only have a uniquely aligned central child at every subsequent level, identified by adding $0$ as a suffix to the code at the current level.

Figure 5. Hierarchical partitioning and encoding of the combined structure. (a) level 0, (b) level 1, and (c) Z curve of the Morton code at level 2. The figures are visualized in Wolfram Mathematica 13.2.1.

Therefore, hierarchical encoding rules are designed, as shown in Figure 6. The encoding form is easy to understand, referred to as the Rhombic Triacontahedron Hexagonal Quadtree on Combined Structure (RTHQCS). The grid level is referred to as $n$, and one Morton digit is added as a suffix to the code when the level increases by one. The first digit ($S$) identifies the number of combined structures, $S=0,1\cdots 9$. The second digit ($R$) identifies the type of initial cells: $0$ represents the hexagonal cell on the upper rhombic surface; $1$ represents the pentagonal cell on the middle rhombic surface; $2$ represents the hexagonal cell on the lower rhombic surface; $3$ represents the pentagonal cell at the north and south vertices of RT. The third to $n+2$th Morton digits, $M[r](r=1,2\cdots ,n)$, identify the quadtree children at each level: $0$ represents the central child; $1$ and $2$ represent the boundary child; $3$ represents the boundary child or external child.

Figure 6. Encoding rules for RTHQCS.

Figure 5(c) shows the Z curve of the Morton code at the second-level grid; thus, the code of each cell is unique, and each digit of the code is related to the geographical location. Assuming that Figure 5 is a diagram of the No. 0 combined structure, the leftmost gray cell in Figure 5(c) is located on the middle rhombic surface, then the Morton digits for the two consecutive partitions are both $3$. Therefore, this cell is coded as $0133$.

This scheme is an advancing quadtree partitioning and Z-order encoding, which enables seamless and non-overlapping grid partitioning globally with identical cell arrangements within 10 combined structures. The pentagonal cells at the south and north vertices are marked in red. Figure 7 shows the results of the grid generation at the second and third levels.

Figure 7. Results of the global grid generation. (a) level 2 and (b) level 3. The figures are visualized in Wolfram Mathematica 13.2.1.

Figure 8. Flowchart showing the transformation from the three-axis integer coordinates to Morton code.

The proposed scheme has only 12 pentagonal cells at each level, which are located at fixed vertices of the combined structure. Their area is five-sixth of the hexagonal cells at the same level. Within each combined structure, the number of cells is identical and there is only one pentagonal cell. The number of hexagonal cells, $C_n$, of the grid system can be calculated as follows: (1) $C_n=10\times (3\times 2^n\times 2^n-1)=30\times 4^n-10$(1) The sphere radius is defined as $R=6371.007180918475\mathrm{km}$, with the same surface area as that of the WGS84 ellipsoid. Table 1 lists the total number of cells at each level and the geometric parametric statistics of the hexagonal cells. Therefore, when the level reaches 23, the average spacing between the cell centers is as low as $0.529m$, satisfying the requirements of high-precision analysis.

Table 1. Geometric parametric statistics of the cells.

Level	Total number of cells	Average area of cells ($k{m}^2$)	Average spacing between cell centers ($\mathrm{km}$)
0	32	17,002,187.390803	4,434.617725
1	122	4,250,546.847701	2,217.308863
2	482	1,062,636.711925	1,108.654431
3	1,922	265,659.177981	554.327216
4	7,682	66,414.794495	277.163608
5	30,722	16,603.698624	138.581804
6	122,882	4,150.924656	69.290902
7	491,522	1,037.731164	34.645451
8	1,966,082	259.432791	17.322725
9	7,864,322	64.858198	8.661363
10	31,457,282	16.214549	4.330681
11	125,829,122	4.053637	2.165341
12	503,316,482	1.013409	1.082670
13	2,013,265,922	0.253352	0.541335
14	8,053,063,682	0.063338	0.270668
15	32,212,254,722	0.015835	0.135334
16	128,849,018,882	0.003959	0.067667
17	515,396,075,522	0.000990	0.033833
18	2,061,584,302,082	0.000247	0.016917
19	8,246,337,208,322	0.000062	0.008458
20	32,985,348,833,282	0.000015	0.004229
21	131,941,395,333,122	0.000004	0.002115
22	527,765,581,332,482	0.000001	0.001057
23	2,111,062,325,329,922	0.0000002	0.000529

3.3. Transformation between the code and three-axis integer coordinates

Hierarchical code directly records hierarchical information of cells, with the highest efficiency of hierarchical indexing, but its neighbor indexing is relatively slow, which generally requires the use of code operations. In contrast, the advantage of the single resolution integer coordinate is that its neighbor indexing is highly efficient. Although the quadtree Morton code is not convenient for neighbor indexing, it is closely related to the integer coordinates of the cells on the rhombic surface (Bai 2011). In this study, a three-axis integer coordinate system is established for convenient neighbor indexing and cross-surface operations, which is also connected to the projection coordinates of the rhombic surface (Section 3.5). Therefore, it is required to design a fast transformation between the code and three-axis integer coordinates, which is a crucial step in this study, tightly linking all components of the grid system.

3.3.1. Transforming the three-axis integer coordinates to code

The process of transforming the three-axis integer coordinates to code is divided into two steps: transforming to the transitional coordinates and addition of the binary coordinates (Figure 8).

For example, Figure 9 shows the transformation process at the second level, explained in detail by the cells on the upper rhombic surface of the combined structure (Figure 9(a)).

Figure 9. Diagram of transforming the three-axis integer coordinates to Morton code at the second level. (a) cells on the upper rhombic surface, (b) cells on the middle rhombic surface, and (c) cells on the lower rhombic surface.

In the first step, the three-axis integer coordinates are transformed to transitional coordinates for easier transformation to the code. In this regard, first, the origin of the three-axis integer coordinate system was translated to the center of the cell, with all Morton digits at $0$. For instance, as shown in Figure 9(a), this translates the origin of the IJK coordinate system to the center of the cell $(0,0,4)$, which corresponds to the Morton digits, $M=00$. Second, a transitional coordinate system, $\mathrm{ij}$, is established such that the coordinate axis coincides with the edge of the rhombus. The transformation method is shown in Eq. (2(2) $\{\begin{array}{c}i=I\\ j=I-K+\mathrm{num}\end{array}$(2) ), in which $\mathrm{num}=2^n$ represents the number of cells at the current level on one edge of the rhombus, such as $\mathrm{num}=4$ at the second level: (2) $\{\begin{array}{c}i=I\\ j=I-K+\mathrm{num}\end{array}$(2) In the second step, the binary transitional coordinates are added to obtain the quaternary Morton digits. Bai, Zhao, and Chen (2005) established the transformation between the Morton code and row–column coordinates of the quadrangular cells. Inspired by the above study, this study realizes the fast transformation between the transitional coordinate and quaternary Morton digits of hexagonal cells. The binary numbers, $B_i$ and $B_j$, correspond to transitional integer coordinates, $i$ and $j$, respectively. Morton digits, $M$, are obtained by adding binary coordinate components, as follows: (3) $M=B_i+2B_j$(3)

The proof of Eq. (3(3) $M=B_i+2B_j$(3) ) is as follows:

Taking a cell on the $i$ axis as an example, its Morton digits are $M_i$; each digit of $M_i$ is $M_i[r](r=1,2,\cdots n)$. As binary numbers consist of $0$ and $1$, the range of decimal numbers represented by $n$ digit binary numbers is $0\sim 1+2+2^2+\cdots +2^{n-1}=2^n-1$. The range of values for the $i$ coordinate component at level $n$ is also $0\sim 2^n-1$, such that it can be precisely represented by $n$ digit binary numbers. The binary number corresponding to the $i$ coordinate of any cell is $B_i$ and each digit is $B_i[r](r=1,2,\cdots ,n)$. When $B_i[r]$ are all equal to $1$, $B_i$ obtains the maximum value, and cells whose $B_i[0]=0$ and $B_i[0]=1$ account for half, respectively. At that time, if $B_i[0]=0$, $M_i[0]=0$; if $B_i[0]=1$, $M_i[0]=1$. Similarly, no matter $M_i[0]=0$ or $M_i[0]=1$, cells whose $B_i[1]=0$ and $B_i[1]=1$ account for half, respectively, such that the value of $M_i[1]$ can also be determined based on $B_i[1]$. By analogy level by level, the binary numbers of $i$ coordinate components exactly correspond to the special structure of the quadtree. The evaluation of each digit in the Morton code can be simplified, as shown in Figure 10, which satisfies Eq. (4(4) $M_i[r]=B_i[r]$(4) ): (4) $M_i[r]=B_i[r]$(4) For cells on the $j$ axis, the Morton digits are $M_j$. The binary numbers corresponding to the $j$ coordinates of any cell are $B_j$. As shown in Figure 10, if $B_j[r]=0$, $M_j[r]=0$; if $B_j[r]=1$, $M_j[r]=2$. Therefore: (5) $M_j[r]=2B_j[r]$(5)

Figure 10. Simplified correspondence between the binary transitional coordinates and Morton code.

For any cell on a rhombic surface, the Morton digits are $M$, satisfying the parallelogram addition criterion, as shown in Eq. (6(6) $M[r]=M_i[r]+M_j[r]$(6) ): (6) $M[r]=M_i[r]+M_j[r]$(6) Therefore, in the transitional coordinate system, according to Eqs. (4(4) $M_i[r]=B_i[r]$(4) )–(6(6) $M[r]=M_i[r]+M_j[r]$(6) ), the Morton code of the cell satisfies Eq. (7(7) $M=B_i+2B_j$(7) ): (7) $M=B_i+2B_j$(7)

Similarly, according to Figure 9, the formula for transforming the three-axis integer coordinates of the middle and lower rhombic surface cells to Morton code can be derived, as listed in Table 2.

Figure 11. Flowchart showing the transformation from the Morton code to three-axis integer coordinates.

Table 2. Transforming the three-axis integer coordinates to Morton code.

Location of cells	Transforming the three-axis coordinates to transitional coordinates	Binary coordinate addition operation
Upper rhombic surface	$\{\begin{array}{c}i=I\\ j=I-K+\mathrm{num}\end{array}$	$M=B_i+2B_j$
Middle rhombic surface	$\{\begin{array}{c}i=K\\ j=I-K\end{array}$
Lower rhombic surface	$\{\begin{array}{c}i=\mathrm{num}-J\\ j=I-J\end{array}$

3.3.2 Transforming the code to three-axis integer coordinates

As each Morton digit is related to the spatial position of cells, it is evaluated level by level to achieve multiplication and translation of the transitional coordinates. Then, an inverse transformation is performed to obtain the three-axis integer coordinates (Figure 11).

Considering the cells on the upper rhombic surface as an example, the process of transforming the quadtree Morton code to three-axis integer coordinates can be deduced as follows:

1. Starting from the initial cell of level 0, i.e. the black-dashed cell in Figure 12, the transition coordinates are $i=j=0$.

2. For $M[1]\sim M[n]$, each Morton digit is iterated through the following operations to obtain the transitional coordinates. First, the transitional coordinate of the central child at the next level is calculated. The coordinate of the central child is twice that of the current cell (Bai 2011), i.e. $i=2i$ and $j=2j$. Second, if $M[r]=0(r=1,2,\cdots ,n)$, this is a central child, and $i=i$, $j=j$; if $M[r]=1$, this is an upper-left boundary child, and $i=i+1$, $j=j$; if $M[r]=2$, this is a lower-boundary child, and $i=i$, $j=j+1$; and if $M[r]=3$, this is a lower-left boundary child, and $i=i+1$, $j=j+1$.

Figure 12. Diagram of transforming the code to transitional coordinates at the second level.

According to Figure 9, for any rhombic surface in the combined structure, the correspondence between the Morton code and operation of the transitional coordinate and transitional integer coordinate, $T_r(r=1,2\cdots n)$, are listed in Table 3. Therefore, the relationship between the transitional coordinates and the Morton code can be summarized as follows: (8) $(i,j)=2\times (2\times (\cdots (2T_1+T_2)+\cdots )+T_{n-1})+T_n=2^{n-1}{T}_1+2^{n-2}{T}_2+\cdots +2T_{n-1}+T_n$(8) For example, the second-level cell in Figure 12, whose Morton code is $23$. The boundary child at the first level is obtained with a Morton code of $2$, i.e. $i=0$ and $j=1$; the coordinate of the central child of cell $(0,1)$ is $(0,2)$. The boundary child of the next level is then obtained whose $i=1$ and $j=3$ because $M[2]=3$.

(3)	The process of transforming the transitional coordinates to three-axis integer coordinates is an inverse process of that shown in Table 2. The conclusions are presented in Table 4.

Table 3. Correspondence between the Morton code and transitional coordinate.

$M[r]$	Central child?	Operation of the transitional coordinate	$T_r$
$0$	√	None	$(0,0)$
$1$	×	$i+1$	$(1,0)$
$2$	×	$j+1$	$(0,1)$
$3$	×	$i+1,j+1$	$(1,1)$

Table 4. Transforming the transitional coordinates to three-axis integer coordinates.

Location of cells	Transforming equation
Upper rhombic surface	$\{\begin{array}{c}I=i\\ K=i-j+\mathrm{num}\end{array}$
Middle rhombic surface	$\{\begin{array}{c}I=i+j\\ K=i\end{array}$
Lower rhombic surface	$\{\begin{array}{c}I=\mathrm{num}-i+j\\ J=\mathrm{num}-i\end{array}$

3.4. Code operations

3.4.1. Hierarchical code operations

Hierarchical code operations are divided into searching for parents’ codes and searching for children’s codes. According to the law of hierarchical partitioning, when searching for parents, only the last digit of the code should be removed; when searching for children, the digit set $\{0,1,2,3\}$ should be added at the end of the code based on the spatial position between the children and the current cell. Owing to the concise code format and clear affiliation between parents and children, the hierarchical code operations are easy to implement.

3.4.2. Neighbor code operations

The neighbor operation efficiency of hierarchical code is relatively low. Therefore, the code is transformed to three-axis integer coordinates, and then a simple and efficient neighborhood query is performed. However, unlike triangles and quadrangles, the hexagonal cells cross the surface boundary of the polyhedron owing to the special shape of hexagons and cannot be completely located within one surface. Therefore, the cross-surface operation of the hexagonal cells is examined in this study.

As shown in Figure 13, the neighborhood query of the hexagonal cells can be classified into three situations: the hexagonal cell and its neighborhoods are located inside one combined structure (referred to as interior cells); the hexagonal cell is located on the boundary of the combined structure (referred to as boundary cells); and the hexagonal cell is located inside the combined structure, but some of its neighborhoods belong to the adjacent combined structure (referred to as semi-boundary cells).

Figure 13. Diagram of the neighborhood query at the second level.

For interior cells (gray cells in Figure 13), this study does not consider the cross-surface problem, but simply add or subtract 1 from the three-axis integer coordinates and finally transform the three-axis integer coordinates of six adjacent cells to codes. As all the negative integer coordinates are set as 0 in this study, the methods for the neighbor operation within one combined structure are not completely identical and can be divided into the following three categories:

1. When the cell is located between the $I$ and $K$ axes, its coordinates are $(i,0,k)$. The $J$ coordinates of the neighborhoods are all $0$, and the neighborhood query rule is shown in Figure 14(a).

2. When the cell is located on the $I$ axis, its coordinates are $(i,0,0)$; then, the neighborhood query rule is shown in Figure 14(b).

3. When the cell is located between the $I$ and $J$ axes, its coordinates are $(i,j,0)$. The $K$ coordinates of the neighborhood are all $0$, and the neighborhood query rule is shown in Figure 14(c).

For the boundary cells (green cells in Figure 13), neighborhoods are firstly searched within the combined structure using the rules in Figure 14. However, other neighborhoods are located on adjacent combined structures. Therefore, this study examined methods for crossing-structure neighborhood queries. The combined structures numbered 0–4 are classified into the first category and denoted as $S1$; the combined structures numbered 5–9 are classified into the second category and denoted as $S2$, and the rhombic surfaces of type 0, 1, and 2 are represented using (0), (1), and (2), respectively. The neighborhood query of the boundary cells locate at the origin of the three-axis coordinate system is easily accomplished; the other boundary cells are located on $S1(0)$, $S1(2)$, $S2(0)$, and $S2(2)$. There is a certain regularity in the cross-structure neighbor code operations, specifically referring to the replacement law of digits between the codes of the boundary cells within the combined structure and the codes of the semi-boundary cells within the adjacent combined structure. As listed in Table 5, there are three methods for replacing Morton digits: the Morton digits of the original cell are $M[r](r=0,1,\cdots ,n)$ and the Morton digits of the cell after cross-structure conversion are $M{}^{\mathrm{\prime }}[r]$. Considering the cell $(1,0,5)$ in the No. 0 combined structure in Figure 13 as an example, this cell is located in $S1(0)$, with the code $0001$ and Morton code $01$. The four neighborhoods are searched within the combined structure, and the code of neighborhood $(0,0,4)$ is $0000$. Then replace $0$ with $3$ in their Morton codes through the conversion using formula 1 in Table 5. The other two neighborhoods in the No. 1 combined structure are obtained, whose codes are $1031$ and $1033$.

Table 5. Morton digit replacement for the cross-structure neighborhood query.

Number	Location of original cells	Location of cells after conversion	Replacement of Morton digit
1	$S1(0)$	$S1(0)$	$M^{\mathrm{\prime }}[r]=\{\begin{array}{c}3,M[r]=0\\ M[r],\mathrm{else}\end{array}$
2	$S1(2)$	$S2(0)$
3	$S2(0)$	$S1(2)$
4	$S2(2)$	$S2(2)$
5	$S1(0)$	$S1(1)$	$M{}^{\mathrm{\prime }}[r]=M[r]+1$
6	$S1(2)$	$S2(1)$
7	$S2(0)$	$S1(1)$	$M{}^{\mathrm{\prime }}[r]=M[r]+2$
8	$S2(2)$	$S2(1)$

Figure 14. Neighborhood query rules for the interior cells.

For the semi-boundary cells (blue cells in Figure 13), some neighborhoods must consider the cross-structure problem. Similarly, a neighborhood query is performed following the method shown in Figure 14. If the neighborhoods exceed the range of the integer coordinate system, the actual neighborhoods in the adjacent combined structure can be found through the inverse replacement of Morton digits listed in Table 5 and can be implemented easily (not discussed here). Considering the cell $(5,0,3)$ in the No. 1 combined structure in Figure 13 as an example, with the code $1131$. The four neighborhoods are searched within the combined structure, and the code of neighborhood $(4,0,3)$ is $1113$. Therefore, the other two neighborhoods $0002$ and $0020$ in the No. 0 combined structure can be obtained through the inverse conversion of formula 5 in Table 5.

3.5. Transformation between the code and geographic coordinates

Although the geospatial data formats are diverse, they primarily rely on the geographic coordinates. Therefore, for application of DGGS to the organization and management of geospatial data, a mutual transformation process must be established between the code and geographic coordinates. This study uses an upper rhombic surface as an example to elaborate on the transformation. The transformation methods for the other two rhombic surfaces in the combined structure are similar.

3.5.1. Transforming the geographic coordinates to code

The process of transforming the geographic coordinates to code is as follows. First, the geographic coordinates are projected on the rhombic surface of the RT to obtain the projection coordinates $(x,y)$. The projection coordinates are then transformed to the three-axis integer coordinate $(I,J,K)$ based on the geometric characteristics of grid partitioning. Finally, $(I,J,K)$ is transformed to the code using the method described in Section 3.3.1.

In the second step, the process of transformation from $(x,y)$ to $(I,J,K)$ is as follows:

1. Given the projection coordinate $(x,y)$ of any point $P$ on the upper rhombic surface, as shown in Figure 15, the straight-line equation of the $K$ axis is $y=\tan ({58.285}^{\circ })x-\cos ({31.715}^{\circ })L$ in the projection coordinate system. The straight-line equation of the I axis is $y=-\cos ({31.715}^{\circ })L$, where $L$ is the length of the edge of the rhombic surface. The Euclidean distances, $d_I$, between point $P$ and the $I$ axis, $d_K$ between point $P$ and the $K$ axis are calculated. Point $P$ is then decomposed into float coordinate $l_I$ on the $I$ axis and $l_K$ on the $K$ axis through geometric relationships. Here, $\overline{C{C}_1}$ is the distance between the centers of adjacent hexagonal cells along the $K$ and $J$ axes directions and $\overline{C{C}_2}$ is the distance between the centers of adjacent hexagonal cells along the $I$ axis direction: (9) $\{\begin{array}{c}{l}_I=d_I/(\cos ({31.715}^{\circ })\times \overline{C{C}_2})\\ l_K=d_K/(\cos ({31.715}^{\circ })\times \overline{C{C}_1})\end{array}$(9)

2. Here, $m_I$ and $m_K$ are the results of taking integers for $l_I$ and $l_K$; the decimal component of point $P$ on the $I$ axis is $r_I=l_I-m_I$ and $r_K=l_K-m_K$ on the $K$ axis. The values of $r_K$ corresponding to points $A$, $B$, and $C$ are ${\displaystyle \frac{1}{3},{\displaystyle \frac{1}{2},{\displaystyle \frac{2}{3}}}}$, respectively. Based on this, the three-axis integer coordinate can be calculated using Eq. (10(10) $\begin{array}{rl}& \{\begin{array}{c}I=m_I,K=m_K,r_I<\sin ({31.715}^{\circ })r_K+{\displaystyle \frac{1}{2}}\\ I=m_I+1,K=m_K,\mathrm{else}\end{array}\left(r_K<{\displaystyle \frac{1}{3}}\right)\\ & \{\begin{array}{c}I=m_I,K=m_K,r_I<{\displaystyle \frac{1}{2}+R_1}\\ I=m_I+1,K=m_K,r_I>1+R_2\\ I=m_I+1,K=m_K+1,\mathrm{else}\end{array}\left({\displaystyle \frac{1}{3}⩽r_K<{\displaystyle \frac{1}{2}}}\right)\\ & \{\begin{array}{c}I=m_I+1,K=m_K+1,r_I>{\displaystyle \frac{1}{2}+R_1}\\ I=m_I,K=m_K+1,r_I<R_2\\ I=m_I,K=m_K,\mathrm{else}\end{array}\left({\displaystyle \frac{1}{2}⩽r_K<{\displaystyle \frac{2}{3}}}\right)\\ & \{\begin{array}{c}I=m_I+1,K=m_K+1,r_I>{\displaystyle \frac{1}{2}-\sin ({31.715}^{\circ })(1-r_K)}\\ I=m_I,K=m_K+1,\mathrm{else}\end{array}\left({\displaystyle \frac{2}{3}⩽r_K<1}\right)\end{array}$(10) ) where $R_1=\sin ({29.943}^{\circ })\times \left({\displaystyle \frac{1}{2}-r_K}\right)/\sin ({28.342}^{\circ })$ and $R_2=\sin ({93.373}^{\circ })\times \left(r_K-{\displaystyle \frac{1}{2}}\right)/\sin ({28.342}^{\circ })$ are obtained according to the law of sines: (10) $\begin{array}{rl}& \{\begin{array}{c}I=m_I,K=m_K,r_I<\sin ({31.715}^{\circ })r_K+{\displaystyle \frac{1}{2}}\\ I=m_I+1,K=m_K,\mathrm{else}\end{array}\left(r_K<{\displaystyle \frac{1}{3}}\right)\\ & \{\begin{array}{c}I=m_I,K=m_K,r_I<{\displaystyle \frac{1}{2}+R_1}\\ I=m_I+1,K=m_K,r_I>1+R_2\\ I=m_I+1,K=m_K+1,\mathrm{else}\end{array}\left({\displaystyle \frac{1}{3}⩽r_K<{\displaystyle \frac{1}{2}}}\right)\\ & \{\begin{array}{c}I=m_I+1,K=m_K+1,r_I>{\displaystyle \frac{1}{2}+R_1}\\ I=m_I,K=m_K+1,r_I<R_2\\ I=m_I,K=m_K,\mathrm{else}\end{array}\left({\displaystyle \frac{1}{2}⩽r_K<{\displaystyle \frac{2}{3}}}\right)\\ & \{\begin{array}{c}I=m_I+1,K=m_K+1,r_I>{\displaystyle \frac{1}{2}-\sin ({31.715}^{\circ })(1-r_K)}\\ I=m_I,K=m_K+1,\mathrm{else}\end{array}\left({\displaystyle \frac{2}{3}⩽r_K<1}\right)\end{array}$(10)

Figure 15. Diagram of transforming the projection coordinates to three-axis integer coordinates at the second level.

However, for boundary cells whose $\mathrm{IJK}$ coordinates exceed the range, similar to the red-dashed cells in Figure 15, these cells actually belong to the adjacent combined structure and require cross-structure processing. First, $1$ is subtracted from their $I$ coordinate, which is equivalent to translating to an adjacent semi-boundary cell. Then, the $\mathrm{IJK}$ coordinates of this neighborhood are converted to code. Finally, the codes of the boundary cells are obtained based on the inverse conversion of the Morton digit replacement listed in Table 5.

3.5.2 Transforming the code to geographic coordinates

The processes of transforming the code to geographic coordinates are as follows. Using the method in section 3.3.2, the code is transformed to three-axis integer coordinates $(I,J,K)$. The projection coordinate $(x,y)$ of the rhombic surface is then calculated based on $(I,J,K)$, and the geographic coordinate $(B,L)$ of the central cell is obtained through inverse projection.

As shown in Figure 16, the three-axis integer coordinate of any cell $P$ on the upper rhombic surface is $(I,0,K)$. Then, the integer coordinate of central cell $O$ on the upper rhombic surface is $\left({\displaystyle \frac{\mathrm{num}}{2},0,\mathrm{num}}\right)$; the cell center of $O$ is the origin of the projection coordinate system. The steps for converting $(I,J,K)$ to $(x,y)$ are as follows:

1. Calculate the integer coordinates difference between point $P$ and central cell $O$ using Eq. (11(11) $\{\begin{array}{c}\mathrm{dI}=I-{\displaystyle \frac{\mathrm{num}}{2}}\\ \mathrm{dK}=K-\mathrm{num}\end{array}$(11) ): (11) $\{\begin{array}{c}\mathrm{dI}=I-{\displaystyle \frac{\mathrm{num}}{2}}\\ \mathrm{dK}=K-\mathrm{num}\end{array}$(11)

2. Calculate the projection coordinates of point $P$ through geometric relationships using Eq. (12(12) $\{\begin{array}{c}x=\mathrm{dK}\times \overline{C{C}_1}\times \sin ({31.715}^{\circ })-\mathrm{dI}\times \overline{C{C}_2}\\ y=\mathrm{dK}\times \overline{C{C}_1}\times \cos ({31.715}^{\circ })\end{array}$(12) ): (12) $\{\begin{array}{c}x=\mathrm{dK}\times \overline{C{C}_1}\times \sin ({31.715}^{\circ })-\mathrm{dI}\times \overline{C{C}_2}\\ y=\mathrm{dK}\times \overline{C{C}_1}\times \cos ({31.715}^{\circ })\end{array}$(12)

Figure 16. Diagram of calculating the projection coordinates according to three-axis integer coordinates at the second level.

4. Results and discussion

As there is no hierarchical encoding scheme for RT DGGS, this study selects HLQT and HHUT as the comparison objects, which are efficient encoding schemes among the icosahedral aperture 4 hexagonal grid systems. They are used to traverse the globe at 0.2° latitude and longitude intervals, with a total of 1.62 million points. This study tests the average time of the three operations at a similar resolution: transforming the geographic coordinates to code, transforming the code to geographic coordinates, and neighbor code operations. Additionally, the proposed scheme of RTHQCS at level $n$ has the same number of cells as the $n+1$ level HLQT, while it has three times the number of cells as the $n$ level HHUT and less than the number of cells as the $n+1$ level HHUT. To fully validate the advantages of RTHQCS, this study ensures that the grid resolution of RTHQCS is similar to other schemes; the number of RTHQCS cells is not less than the other schemes. Therefore, this study selects the high-level grid of 13–20, ensuring that HLQT is always one level higher than RTHQCS, and HHUT is always the same level as RTHQCS.

Different DGGSs adopt different projection methods; for example, HLQT and HHUT adopt the icosahedral Snyder projection, whereas RTHQCS employs the RT Snyder projection (Liang et al. 2022). In the process of transforming the geographic coordinates to code, the impacts of different projection methods on the transforming efficiency vary. Both the forward and inverse projections occupy a portion of the time. The inverse Snyder projection requires an iterative solution which accounts for the majority of the transformation time from the code to geographic coordinates. Therefore, to fairly compare various encoding schemes, this study eliminates the time required for forward and inverse projections.

The time unit is microseconds. Each operation is tested 10 times, and the average value is used. The efficiency ratio is the ratio of the time consumed in HLQT or HHUT to the time consumed in RTHQCS. All programs are compiled into release versions and tested on a compatible machine (hardware configuration: Inter Core i5-10400F CPU2.90 GHz, 16G RAM. KIOXIA-EXCERIA 480G SSD; operation system: Windows 10 x64 Enterprise LTSC; development tools: Microsoft Visual C++ Enterprise 2022). Tables 6–8 list the results. Figures 17–19 show efficiency comparisons.

Figure 17. Comparison of transforming the geographic coordinates to code. (a) time comparison and (b) efficiency ratio.

Figure 18. Comparison of transforming the code to geographic coordinates. (a) time comparison and (b) efficiency ratio.

Figure 19. Comparison of the neighbor code operation. (a) time comparison and (b) efficiency ratio.

Table 6. Experimental results of transforming the geographic coordinates to code.

Level	RTHQCS($\mathrm{\mu s}$)	HLQT($\mathrm{\mu s}$)	Efficiency ratio	HHUT($\mathrm{\mu s}$)	Efficiency ratio
13	0.20690	4.47440	21.62602	0.74024	3.57777
14	0.20820	4.72862	22.71224	0.77747	3.73428
15	0.21323	5.08861	23.86484	0.83522	3.91707
16	0.21757	5.42810	24.94920	0.87935	4.04174
17	0.22611	5.77383	25.53506	0.93629	4.14077
18	0.23450	6.14448	26.20257	0.98701	4.20901
19	0.24980	6.55426	26.23781	1.05462	4.22181
20	0.25813	6.97257	27.01175	1.10758	4.29076

Table 7. Experimental results of transforming the code to geographic coordinates.

Level	RTHQCS($\mathrm{\mu s}$)	HLQT($\mathrm{\mu s}$)	Efficiency ratio	HHUT($\mathrm{\mu s}$)	Efficiency ratio
13	0.05138	1.90158	37.00945	0.07413	1.44282
14	0.05404	2.09127	38.70003	0.08797	1.62788
15	0.05762	2.22107	38.54757	0.09880	1.71478
16	0.05969	2.33763	39.16473	0.10458	1.75215
17	0.06262	2.49769	39.88462	0.11824	1.88806
18	0.06482	2.60020	40.11546	0.12640	1.95010
19	0.06303	2.74590	43.56770	0.13297	2.10972
20	0.06171	2.89198	46.86472	0.14592	2.36470

Table 8. Experimental results of the neighbor code operation.

Level	RTHQCS($\mathrm{\mu s}$)	HLQT($\mathrm{\mu s}$)	Efficiency ratio	HHUT($\mathrm{\mu s}$)	Efficiency ratio
13	0.26854	8.46650	31.52755	1.04953	3.90825
14	0.28199	9.33151	33.09150	1.12717	3.99718
15	0.29548	10.35248	35.03615	1.21314	4.10564
16	0.31648	11.20055	35.39069	1.30505	4.12359
17	0.32935	12.21800	37.09786	1.38279	4.19860
18	0.34931	13.27694	38.00893	1.48737	4.25800
19	0.36780	14.58361	39.65114	1.58159	4.30015
20	0.38438	15.74993	40.97521	1.67454	4.35650

Experimental results indicate that:

1. RTHQCS has the highest efficiency in transforming the geographic coordinates to code (Figure 17(a)); the efficiency advantage gradually improves with an increase in level (Figure 17(b)), reaching 27.01-fold that of HLQT and 4.29-fold that of HHUT (Table 6). RTHQCS transforms the three-axis integer coordinates to code, and it is a fast algorithm, almost unaffected by the level. Because this transformation is achieved directly through the addition of binary transitional coordinates, while complex judgments related to the level are avoided. In contrast, both HLQT and HHUT require slow code addition operations using individual digits one by one to obtain the code. As the level increases, the length of the code increases, which increases the computational complexity of the code addition. Additionally, HLQT consists of 12 vertex tiles and 20 face tiles, which is a more complex structure and has the lowest efficiency.

2. RTHQCS has the highest efficiency in transforming the code to geographic coordinates (Figure 18(a)); the efficiency advantage gradually improves with an increase in level (Figure 18(b)), reaching 46.86-fold that of HLQT and 2.36-fold that of HHUT (Table 7). RTHQCS code is composed only of $0,1,2,3$; when the digit is $0$, it does not need to participate in the calculation (Table 3). Additionally, the projection coordinates are obtained in the plane using simple geometric relationships after the code is transformed to three-axis integer coordinates. In contrast, the digits of the HLQT code at each level must participate in the calculations; its first digit is complicated, consuming significant processing time. Although the HHUT code does not need to participate in the calculation when the digit is $0$, it can only obtain the tiles’ Cartesian coordinates after scaling and accumulation of complex number codes; it still requires multiple judgments to obtain the projection coordinates.

3. RTHQCS has the highest efficiency in neighbor code operation (Figure 19(a)); the efficiency advantage gradually improves with an increase in level (Figure 19(b)), reaching 40.97-fold that of HLQT and 4.35-fold that of HHUT (Table 8). This study combines three adjacent rhombic surfaces into a combined structure and establishes a three-axis integer coordinate system to simplify the problem of cross-surface neighborhood queries. Adjacent cells within the combined structure can be obtained easily by adding or subtracting 1 in the integer coordinates (Figure 4), while the cross-structure problem between adjacent combined structures can be solved directly through the replacement of code digits (Table 5), which avoids rotation and translation in traditional cross-surface operations. The neighbor code operations of HLQT are achieved through code addition operations with six directional vectors, with the lowest efficiency. In the HHUT scheme, some neighborhoods can be directly obtained by changing the end digit, but other neighborhoods still require code addition operations. Both the comparison schemes require searching the code addition operation lookup table; as the level increases, the length of the code increases and the code digits are processed one by one in code addition, which reduces the computational efficiency. In conclusion, the neighbor queries based on integer coordinate is much quicker than code addition and unrelated to grid level. Therefore, RTHQCS has the highest efficiency in the neighbor code operation and is not affected by the level.

In summary, compared to HLQT and HHUT, RTHQCS not only enables the construction of spherical grids with superior geometric properties but also has significant efficiency advantages in the transformation between the code and geographic coordinates as well as neighbor code operations.

5. Conclusion

Based on the current situation for RT DGGS, multi-resolution hierarchical encoding schemes are lacking and frequent cross-surface operations reduce the efficiency of data processing. Considering the structural characteristics of the RT aperture 4 hexagonal grids, this study proposes a hierarchical partitioning and encoding scheme, designs code transformation and code operation algorithms, and verifies the superiority of the proposed scheme through comparative experiments.

The proposed RTHQCS has broad application prospects. First, the method can obtain global grids with smaller distortions, offering a higher degree Earth fitting. Therefore, the geographic locations of the grid cells are closer to the real world, beneficial for accurate data modeling. Second, the fast transformation between the code and geographic coordinates enables RTHQCS to organize and analyze geospatial data efficiently, such as the integration of terrain data and environmental modeling. Third, the fast neighbor code operations provide potential for efficient spatial analysis, such as path planning.

This study only designs the nearest neighbor operation, which can obtain the nearest 6 hexagonal cells of a specific cell. However, the nearest neighbor operation is insufficient for offsetting regions, such as in the case of looking for sharing-bike within 500 meters. In the future work, it is supposed to develop a ringed-region neighbor operation algorithm to accelerate the process of searching for cells within a ring.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available in figshare.com with the identifier(s) https://figshare.com/s/7df09448a64226b592ec.

Additional information

Funding

This work was supported by the National Defense Science Innovation Special Zone Program of China [grant number: 20-163-14-LZ-001-003-01] and the Special Science Fund for Innovation Ecosystem Construction of National Supercomputing Center in Zhengzhou [grant number: 201400210100].