The Transform between Geographic Coordinates and Location Codes of Aperture 4 Hexagonal Grid System

Discrete global grid system is a new data model which supports the fusion processing of multi-source geospatial data. In discrete global grid systems, all cell operations can be completed by codes theoretically, but most of current spatial data are in the forms of geographic coordinates and projected coordinates. It is necessary to study the transform between geographic coordinates and grid codes, which will support data entering and getting out of the systems. This paper chooses the icosahedral hexagonal discrete global system as a base, and builds the mapping relationships between the sphere and the icosahedron. Then an encoding scheme of planar aperture 4 hexagonal grid system is designed and applied to the icosahedron. Basing on this, a new algorithm of transforms between geographic coordinates and grid codes is designed. Finally, experiments test the accuracy and efficiency of this algorithm. The efficiency of code addition of HLQT is about 5 times the efficiency of code addition of HQBS.


INTRODUCTION
Discrete Global Grid Systems (DGGS) is a new global-oriented data model, of which the main characteristic is the digital representation of geospatial location.DGGS partition the globe into uniform grid cells in hierarchical structure and address every cell using location code to substitute for traditional geographic coordinates in operation.Compared with the traditional spatial data models, DGGS take the globe as the research object and provide worldwide uniform space datum.
Each grid cell is correspond to one geographical location and cell size will change in terms of the level of grid, which help process of multi-resource, isomerous and multi-resolution geospatial data.Moreover, in DGGS, geometrical operation of grid cells can be achieved completely by location codes which improve efficiency of data computation and process. 1   Common cell shapes used by DGGS are triangles, quadrangles and hexagons.Compared with other two shapes, the topological relationship of hexagonal grids is consistent, which is beneficial for spatial analysis like adjacency and connection (Sahr et al., 2003;Ben et al., 2015a).Spatial sampling efficiency of hexagonal grids is highest and helps data visualization (Sahr, 2011).The processing efficiency of sampling values is about 20 percent to 50 percent higher than quadrangles (Sahr, 2011).
In the regular polyhedral DGGS, the icosahedron has the most face, hence its geometric distortion after projection is minimum.Now, many researchers at home and abroad have made relevant research into the icosahedral hexagonal discrete grid systems.Sahr et al. (Sahr, 2005;Sahr, 2008) realized the fast indexing of cells at different resolutions based on the Icosahedral Snyder Equal Area aperture 3 Hexagonal grid (ISEA3H).Vince and Zheng (Vince and Zheng, 2009) designed the algorithms of location code operation and Fourier Transform of the icosahedral aperture 3 hexagonal grid system.Tong et al. (Tong , 2010;Tong et al., 2013) designed an indexing structure called the Hexagonal Quad Balances Structure (HQBS).Ben et al. (Ben et al., 2015a;Ben, 2005;Ben et al., 2007;Ben et al., 2010;Ben et al., 2011;Ben et al., 2015b) researched into the generation algorithm, code indexing and real-time display of hexagonal DGGS.In engineering, PYXIS Innovation Corporation in Canada developed the spatial data integration software based on the ISEA3H grid system and its core technology is PYXIS grid indexing patent (Peterson, 2006).
Although current research has made considerable progress, some deficiency exists.In the PYXIS scheme, the direction of grid cells will change with levels of resolution, which is not beneficial to spatial analysis like adjacency.The concept model of HQBS is so complicated that there are frequent failure in code regularization and operaions will rollback, which seriously affects computational efficiency.The operations of cells can be completely realized by location code in grid systems, but at present most spatial data are still represented and stored in the form of traditional geographic coordinates or projection coordinates.To ensure that spatial data can transform between the traditional data organization framework and grid systems at a high speed, it is necessary to make research on conversion between geographic coordinates and location codes.Above all, this paper chooses the icosahedral DGGS and firstly establishes the mapping relationship between the regular icosahedron and the sphere.Then, an encoding scheme of the planar aperture 4 hexagonal grid system is designed and applied into the icosahedron skillfully.Based on this, we desinged a new algorithm of conversion between geographic coordinates and location codes of grids and compares with the scheme HQBS to verify validity and high efficicency of the algorithm.

LOCATION OF THE REGULAR ICOSAHEDRON AND MAPPING RELAITON WITH THE SPHERE
In order to establish the mapping relation between the surface of a regular icosahedron and the sphere, firstly we should determine the location relationship between the regular icosahedron and the earth spheriod.Locate two vertices of the regular icosahedron at the north and south poles respectively and locate another one vertex at (0 ± ;25:56505 ± N ).Then flatten the regular icosahedron and number each vertex or assign every vertxe an index as illustrated in Figure 1.The rule is as the following: Vertex located at the North pole and South pole are indexed as 0 and 11 respectively.The number 1 vertex is at the baisi meridian and the restcan be done in the same manner.Apart from vertices, faces are needed to be indexed as well: faces north of 25:56505 ± N are assigned 0-4 and faces sharing common edges with them are assinged 5-9.Similarly, faces south of 25:56505 ± S assigned 15-19 and faces sharing common edges with them are assinged 10-14.According to the above relations, given a point P (B ;L ) whose geograohic coordinates, that is longitute and latitude, is known, the index T P of the triangular face P (B ;L ) belongs to can be obtained.
For the projection between the regular icosahedron and the sphere, this paper adopts Snyder Equal-area Map Projection for polyhedral globes to establish the mapping relations of single face of the icosahedron and the sphere.The Snyder projection ensures the continuity of network of longitude and latitude and decrease the distortion at the same time.The forward and inverse solution, detailed procedure and algorithms of Snyder Equal-area Map Projection can be referenced in Snyder (1992) and Ben et al. (2006) but not discussed in this paper.
In terms of the above location and projection, given a point P (B ;L ), compute the triangular face number T P it locates at, then via Snyder projection, the 2D Cartesian coordinate

Encoding scheme and code addition of the planar A4H
Above analysis, in order to complete the conversion, we need to research relationship between locaiton codes on the triangular faces and the correspongding 2D Cartesian coordinates.Centers of grid cells are called lattice, which are identified with grid cells in our research.For convnence of code substitution, this paper adopts complex numbers to represent lattice locations (Vince, 2006a(Vince, , 2006b)).Thereinaftet, let n denote the partition level of grids, and the higher the level, the finer the resolution or the size of grid cell.According to geometry of the aperture 4 hexagongal grid system (A4H), lattice at level comprises lattice and midpoints of each cell edges at level n ¡ 1.
On the complex plane, let i and D = f0;!;! 2 , ;! 3 g, the set of lattice of A4H is given by where '
In equation ( 1), representation of each lattice is unique, which satifies the requirement of the location code for uniqueness.Let n are all complex numbers so it is not convenient for expression and computation.Now, replace d i , 1 6 i 6 n , with location codes.Omit the base of power exponents in equation ( 1) and make replacement as below ;2 6 i6 n), where the digit set D 0 becomes E 0 = f0 ;1 ;2 ; 3 ;4 ;5 ;6 g and similarly D becomes E = f0 ;1 ;2 ;3 g.
Hence, we get the unique identifier of  Code operation is more suitable for computaionally processing geometric operation of cells and used to substitute for traditional operation using floating points.In HLQT, the location code records the distance and direction of the lattice to the origin whose code is conpletely made up of digit 0 and the addition and substraction of location codes are equivatent to the addition and substraction of vectors on the complex plane and satisfy parallelogram law of vectors.In the conversion from geographic coordiantes to locaiton codes, we also use the addition operation of location codes, hence we simply introduce some properties of code addition and operation rule by induction.
In HLQT, Code set H n ;n > 1, is a commutative group G fH n ;© g and has the following properties: 1. Closure.
123 © 010 = 233 .And these properties make G fH n ;© g the Abelian group.The operation rule of code addtion is similar to the decimal addition essentailly.In operation , we reference the addition table as Table 1.

Extend of the encoding scheme onto the icosahedron
Hereinbefore, we discussed the encoding scheme on the infinite 2D plane, but the surface of the regular icosahedron is closed.
In order to extend the encoding scheme to the regular icosahedron, 'face tile' and 'vertex tile' are used to depict its , where a, b and c are indexes of the vertices of the triangle and i is the index of the triangle, i = 0;1;¢¢¢ ;19.T a ;b;c i comprises a complete P i and part of V a , V b and V c .According to Section 2, we can compute T P , which is i herein equivalently, using geographic coordinates.If T P is obtained, it is easy to get the indexed of vertices in term of Figure 1 and the code structure in this triangular face.
The arrangement of location codes on the face tile is consistent with that on the vertex tile except the number of grid cells as illustrated in Figure 4, which has no effect on the algorithm.
According to this, make the tile as the unit when executing code operation.Hence define the 2D coordinate system O 0 ¡ X 0 Y 0 in the face tile and vertex tile.As shown in Figure 4, define the lattice whose code is all made up of digit 0 as the origin O 0 , the line of origin and lattice whose code is 010 as X-axis and Y -axis conforms to the right-handed coordinate system.
Obviously, O 0 ¡ X 0 Y 0 of the face tiles is equal to the corresponding O ¡ X Y in the triangular faces, but O 0 ¡ X 0 Y 0 of the vertex tile will translate and rotate according to the triangular face and vertex it associates with.

COORDINATE AND LOCATION CODE
The conversion has two forms.One is from geographic coordinates to location codes, another is from location codes to geographic coordinates.

Conversion from geographic coordinates to location codes
In this conversion, we introduce a three-axis coordinate system IJK ¡ O as shown in Figure 5.The origin of IJK O is consistent with O 0 ¡ X 0 Y 0 , and J-axis is the same as Y 0 -axis.
The included angles of the three axes in the positive direction are all 120 ± , which divides the plane into six quadrants.
The basic idea of this conversion is that given a point P (B ;L ) on the sphere, compute the 2D Cartesian coordinate p(x p ;y p ) on the O ¡ X Y of a single triangular face according to the location parameters of the regular icosahedron and forward solution of Snyder Projection.Then, in terms of the relationship of the triangular faces, face tile and vertex tile, converse p(x p ;y p ) to p(x q ;y q ) on the correspongding tile O 0 ¡ X 0 Y 0 this point belongs to.In IJK ¡ O , using the parallelogram law, project p(x q ;y q ) to the 2 axes wich are more closer to it and obtain two component code ® 1 and ® 2 in the two axes.Finally, add ® 1 and ® 2 in terms of rule of code addition and get the location code where p is located.
The (A i ;® ) or (V i ;® ) comprises the tile and location code associated with p.The choose of projection axes can be judged by the angle between the line of p and the origin O and X 0 axis.I , J and K axes are lines through a series consecutive lattice and via induction, the lattter n ¡ 1 digits of location codes of lattice at these three axes have the following relaitonship with k > < > : where k denotes the position of lattice at the axis and D B in (k) represents the binary number of k.Although the first digit doesn't match the rule of binary system, its relation with k is easy to obtain as well.Finally, in terms of addition, the location code of the cell p is inside is given by

Conversion from location codes to geographic coordinates
The conversion from location codes to geographic coordinates is easier.For the icosahedral code structure, the same code will exist in different tiles, hence apart from the location code, the index of the face tile or vertex tile is needed to be given as well.
Given (A i ;® ) or (V i ;® ), the detailed steps are as the following: 1. Compute the complex number representation of the location code ® .
According to Section 2, the lattice is number of a special form on the complex plane and the form of location code is just another efficient representation for it.Location code and complex number of the lattice can transform each other.Let the level of grid be n , = e 1 e 2 ¢¢¢e n , ² e n = 0 ¢¢¢0 | {z } n ¡ 1 ;n > 2; e 2 f0 ;1 ;2 ;3 g, the complex number of ® denoted C (® ) is that here the real part and imaginary part of C (® ) are corresponding to (x q ;y q ) in the coordinate system of tiles .
2. Translate and rotate (x q ;y q ) into p(x p ;y p ) in the in the coordinate system of triangular faces.

Via the inverse computation of Snyder Equal-area Map
projection, obtain the resulitinf (B ;L ).The detailed process of inver computation can refer to Snyder (1992) and Ben et al. (2006).

EXPERIMENTS AND ANALYSIS
To verify the efficiency of the algorithm, we test the efficiency of code addition for the HLQT scheme and compare it with HQBS in that both the 2 shcemes use A4H and Snyder projection and the important factor influencing efficiency is the effiency of code operation.We adopt Visual C++ 2012 ad the development tool, compile the two algorithms into release versions and operate on the same compatible PC (Configuration: Intel Core i5-4590 CPU@3.3GHz,8G RAM, KingSton 120GB SSD; Operating System：Win 7 X64 Ultimate SP1; ).The result and the correlation curve of efficiency are shown in Table 3 and Figure 6 respectively.We can see from Figure 6 that the efficiency of code addition of HLQT is obviously higher than HQBS and about 5 times the efficiency of code addition of HQBS.
Figure 1.Indexes of vertices and faces in the icosahedron and the coordinate system in triangular face in Figure 2. Lattices at different levels and their relations are expressed using circles and triangles with different size.For convenience of representation, thereinafter name the hierarchical structure of lattice of A4H Hexagonal Lattice Quad Tree (HLQT).

Figure 2 .
Figure 2. Representation of codes in the 3 rd level and examples of code operation structure.The flattened icosahedron is shown in Figure 3.The face tiles are centered at the center of each triangle of the icosahedron and the cells of them are filled in white.The vertex tiles are centered at the center of each vertex of the icosahedron and the cells of them are filled in red and green.Denote the set of lattice of the icosahedral A4H at level n G n (n > 1).G n comprises two parts.One part is 20 face tiles A n ;i = 0;1;¢¢¢ ;19 and i is the index of the triangular faces the face tile is centered at.Another part is 12 vertxe tiles V j n ;j = 0;1;¢¢¢ ;11 and j is the index of the vertex the vertex tile is centered at.That is

Figure 3 .
Figure 3. Vertex tiles and face tiles in the flattened icosahedron Figure 5. Conversion from geographic coordinates and location codes

Figure 6 .
Figure 6.Comparison of computational efficiency about addition in the two schemes