CN117689840A - Coding and operation method and device of rhombic icosahedron global discrete grid system - Google Patents

Abstract

The invention discloses a coding and operation method and a device of a rhombic icosahedron global discrete grid system. And secondly, a set of brand-new rhombic icosahedron hierarchical segmentation and quadtree coding schemes are established based on the combined structure, and rapid interconversion of codes and triaxial integer coordinates is constructed. Then, the coding operation method of the scheme is provided, hierarchical query is easy to realize, the inner unit of the combined structure directly utilizes the three-axis integer coordinates to realize the adjacent query, and the method that the boundary unit directly utilizes the code element to replace to realize the adjacent query is generalized according to the coding characteristics. Finally, the precise conversion between the geographic coordinates and the codes is realized through the strict geometrical relationship between the projection coordinates and the triaxial integer coordinates.

Description

Coding and operation method and device of rhombic icosahedron global discrete grid system

Technical Field

The invention relates to the technical field of geographic information, and can be applied to organization management and fusion analysis of wide-area multi-source heterogeneous geospatial data, in particular to a coding and operation method and device of a rhombic icosahedron global discrete grid system.

Background

Currently, the data processing capability of the traditional geographic information system (Geographic Information System, GIS) is difficult to meet the requirements due to large data characteristics such as wide coverage range, huge data volume, various types and formats, complex data structure and the like of the geographic space data, and a new scheme for supporting global mass geographic space data fusion processing is needed. The global discrete grid system divides the earth space into a seamless and non-stacked multi-scale space hierarchical structure, is similar to a spreadsheet covering the earth space, and is a space-time data unified organization, processing and mining framework with great potential. The traditional theodolite grid has serious deformation of units in high-latitude areas, and a large amount of data redundancy exists when the traditional theodolite grid is used for data organization. The polyhedral global discrete grid system adopts a specific method to recursively split the surface of the polyhedron, and then projects the polyhedron surface to the surface of the earth to form a global continuous and approximately regular uniform spherical hierarchical grid structure, so that the polyhedral global discrete grid system becomes a hot spot for current research and application.

The grid subdivision cell types are generally classified into triangles, quadrilaterals, and hexagons, where the adjacency consistency, highest spatial sampling rate, and angular resolution of the hexagons make them more desirable geometric properties. The split aperture refers to the ratio of the areas of adjacent hierarchical units, and the aperture of the hexagonal grid is 3 of three holes, four holes and seven holes. The three-hole grid can provide more resolution, and the transition between adjacent layers is the smoothest; however, the directions of adjacent hierarchical units are changed periodically by 30 degrees, so that the arrangement rules of the odd-even layer units are different, and the hierarchical coding operation is complex. The seven-hole grid can achieve approximate consistency, and the unit membership problem is well solved; but the area of the cells of the adjacent layers of the seven-cell network is too much changed, and the direction is rotated by 19.1 degrees. The directions of adjacent hierarchical units of the four-hole grid are unchanged, so that a multi-layer structure is convenient to build. White et al statistically analyze the effect of five Berrader polyhedrons on the deformation property of the grid, and find that the deformation caused by the regular icosahedron is the smallest, so most of the current researches are the regular icosahedron hierarchical coding scheme, and the common open source libraries such as DGGRID, openEAGGR, geogrid all adopt the regular icosahedron to construct the DGGS.

Recently, students try to use a grid system with more faces and more excellent geometric properties of a Catalan polyhedron structure, and develop a brand new DGGS design idea. The Hall et al use one hundred twenty faces to construct DGGS, but the single faces are irregular triangular faces, so that the DGGS is only suitable for constructing triangular grids and has complex hierarchical relationship; liang et al construct a rhombic icosahedron hexagonal grid, each surface is a gold diamond with the same size, the diamond is suitable for the subdivision of various units, a group of adjacent sides of the diamond surface are natural row and column coordinate systems, the positions of the units are easy to describe, and the advantages of the diamond icosahedron hexagonal grid system in terms of compactness in shape and grid deformation are proved; wang et al realized rapid grid generation of irregular areas by using rhombic icosahedron SD isostock projected hexagonal grids.

Disclosure of Invention

The existing rhombic icosahedron grid system coding schemes all adopt integer coordinate schemes with single resolution, the deep research on multi-resolution coding and operation thereof is lacking, and continuous multi-level grid analysis cannot be supported; and the cross-plane operation among thirty diamond planes is frequent, which limits the application of the grid system in large-area geospatial data processing and analysis. The invention provides a coding and operation method and device of a rhombic icosahedron global discrete grid system (Discrete Global Grid Systems, DGGS) aiming at the problems, and the obtained global discrete grid system can be applied to geospatial data organization and storage, indexing and calculation. The invention designs a layering subdivision and coding method of hexagonal grid cells on a rhombic icosahedron based on simplified cross-plane operation, and designs operations such as rapid adjacent query and layering query of cell coding, mutual conversion of coding and geographic coordinates.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

the invention provides a coding and operation method of a rhombic icosahedron global discrete grid system, which comprises the following steps:

step 1, dividing a rhombic icosahedron into 10 combined structures, wherein the serial numbers of the combined structures are 0-9, and each combined structure is formed by combining adjacent upper, middle and lower rhombic faces; establishing a triaxial integer coordinate system in each combined structure;

step 2, for each combined structure, constructing an initial unit comprising an upper diamond-shaped hexagonal unit, a middle diamond-shaped pentagonal unit and a lower diamond-shaped hexagonal unit, wherein the number of the upper diamond-shaped hexagonal unit is 0, the number of the middle diamond-shaped pentagonal unit is 1, and the number of the lower diamond-shaped hexagonal unit is 2; in addition, the initial pentagonal units at the south and north vertexes of the rhombic icosahedron are numbered 3 and respectively assigned to a 2-7 combined structure; different four-hole hierarchical splitting and coding methods are adopted for different initial units;

step 3, constructing the interconversion of the hierarchical coding and the triaxial integer coordinates;

step 4, designing a coding operation method;

and 5, constructing the mutual conversion of the geographic coordinates and the grid codes so as to integrate the geographic space data into the grid system and obtain the rhombic icosahedron global discrete grid system.

Further, the four-hole hierarchical subdivision and coding method for different initial units comprises the following steps: for the initial pentagon unit of the north and south vertexes, only a unique sub-unit with aligned centers is identified by a code element 0 after each layering.

Further, the four-hole hierarchical subdivision and coding method for different initial units comprises the following steps:

performing four-hole hierarchical segmentation and coding on other initial units except for the initial pentagon unit of the north-south vertex, representing subunits at different positions by 0-3 code elements, wherein 0 is used for identifying a central subunit, 1-2 is used for identifying a boundary subunit, 3 is used for identifying a boundary subunit or an external subunit, and the hexagonal units of the upper diamond surface and the hexagonal units of the lower diamond surface and the two subunits adopt a hexagonal quadtree segmentation and Morton code coding scheme, so that 3 boundary subunits are adjacent to the current unit; and carrying out hierarchical subdivision coding on the pentagon unit and the sub-units of the pentagon unit with the middle diamond surface, so that 1 sub-unit is completely positioned outside the current unit in the space position.

Further, the step 3 includes:

step 3.1, converting the triaxial integer coordinates into codes, and for the diamond-shaped surface on the combined structure, firstly, carrying out coordinate transformation on the triaxial integer coordinates (I, J, K) of the combined structure according to the formula (1) to obtain diamond-shaped surface transition coordinates (I, J):

wherein num=2 ⁿ Representing the number of units on one side of the diamond-shaped surface when the current grid is layered;

the binary components of the transition coordinates are then added according to equation (2) to obtain Morton codes:

M＝B _i +2B _j (2)

wherein B is _i And B _j Binary numbers corresponding to the transition coordinates i and j of the units are respectively represented;

and 3.2, converting the codes into triaxial integer coordinates, judging Morton codes layer by layer according to the relation between the Morton codes and the unit space positions for the diamond surface on the combined structure, multiplying and translating the transition coordinates, and finally obtaining the triaxial integer coordinates through inverse transformation.

Further, the step 3.2 includes:

step 3.2.1, setting the transition coordinates of the initial grid unit of the 0 th layer to be (i, j) = (0, 0);

step 3.2.2, executing the following operations on each Morton code to finally obtain the transition coordinates of the unit: calculating transition coordinates i=2i, j=2j of a central subunit of the next layer, and if the r-th Morton symbol Mr of the subunit is 0, i.e. the unit is the central subunit, i=i, j=j; if M [ r ] =1, i.e. the cell is the upper left boundary subunit, i=i+1, j=j; if M [ r ] =2, i.e. the cell is the lower boundary subunit, i=i, j=j+1; if M [ r ] =3, i.e. the cell is the lower left boundary subunit, i=i+1, j=j+1;

and 3.2.3, converting the transition coordinates into three-axis integer coordinates, wherein the process is the inverse process of the formula (1).

Further, the step 4 includes:

step 4.1, designing a hierarchical coding operation method, namely searching for a parent unit code and searching for a child unit code, and removing the tail bit of the code when searching for the parent unit code; when searching for the subunit codes, adding a code element set {0,1,2,3} in the tail bits of the codes according to the space positions of the subunit and the father unit;

step 4.2, designing a neighboring encoding operation method, wherein the neighboring operation of the hexagonal unit is divided into three cases: the hexagons and adjacent units thereof are all positioned inside the affiliated combined structure, the hexagons are positioned on the boundary of the affiliated combined structure, the hexagons are positioned inside the affiliated combined structure, but part of adjacent units thereof are affiliated to the adjacent combined structure, the units corresponding to the first case are called as internal units, the units corresponding to the second case are called as boundary units, and the units corresponding to the third case are called as half-boundary units.

Further, the step 4.2 includes:

step 4.2.1, completing the adjacent operation of the internal unit by adding or subtracting 1 to the three-axis integer coordinate;

step 4.2.2, for boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, wherein other adjacent units are positioned on the adjacent combined structure, classifying the combined structure with the number of 0-4 into a first class, marking the combined structure as S1, correspondingly marking the combined structure with the number of 5-9 as S2, respectively using (0), (1) and (2) to represent upper, middle and lower diamond surfaces of the combined structure, and inquiring the adjacent combined structure half boundary units by the boundary units according to different unit positions to obtain 3 Morton code element replacing methods as follows:

wherein M r represents Morton symbols of the unit before transformation, and M' r represents Morton symbols after transformation;

step 4.2.3, for half boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, and then finding other adjacent units on the adjacent combined structure by the reverse process of symbol substitution.

Further, the step 5 includes:

step 5.1, converting the geographic coordinates into codes: firstly, projecting a geographic coordinate to an icosahedron diamond-shaped surface to obtain a projection coordinate (x, y); then converting the projection coordinates into three-axis integer coordinates (I, J, K) according to the grid subdivision geometric characteristics; finally, converting the three-axis integer coordinates into codes by the method of the step 3.1;

step 5.2, converting the codes to geographic coordinates: firstly, converting codes into three-axis integer coordinates (I, J, K) by a method of the step 3.2; then calculating diamond plane projection coordinates (x, y); finally, the geographical coordinates (B, L) of the cell centers are obtained through back projection calculation.

Further, after the step 5, the method further includes:

and 6, performing quality evaluation and comparison on the established rhombic icosahedron global discrete grid system.

Another aspect of the present invention provides a coding and computing device for a rhombic icosahedron global discrete grid system, including:

the first construction module is used for dividing the rhombic icosahedron into 10 combined structures, the serial numbers of the combined structures are 0-9 in sequence, and each combined structure is formed by combining an upper diamond surface, a middle diamond surface and a lower diamond surface which are adjacent; establishing a triaxial integer coordinate system in each combined structure;

the hierarchical splitting and encoding module is used for constructing an initial unit for each combined structure, wherein the initial unit comprises an upper rhombic-surface hexagonal unit, a middle rhombic-surface pentagonal unit and a lower rhombic-surface hexagonal unit, the upper rhombic-surface hexagonal unit is numbered to be 0, the middle rhombic-surface pentagonal unit is numbered to be 1, and the lower rhombic-surface hexagonal unit is numbered to be 2; in addition, the initial pentagonal units at the south and north vertexes of the rhombic icosahedron are numbered 3 and respectively assigned to a 2-7 combined structure; different four-hole hierarchical splitting and coding methods are adopted for different initial units;

the second construction module is used for constructing the interconversion between the hierarchical coding and the triaxial integer coordinates;

the coding operation design module is used for designing a coding operation method;

and the third construction module is used for constructing the mutual conversion between the geographic coordinates and the grid codes so as to integrate the geographic space data into the grid system and obtain the rhombic icosahedron global discrete grid system.

Further, the method further comprises the following steps:

the quality evaluation and comparison module is used for evaluating and comparing the quality of the established rhombic icosahedron global discrete grid system.

Compared with the prior art, the invention has the beneficial effects that:

(1) The rhombic icosahedron has higher fitting degree with the earth, can construct a grid system with smaller deformation and better property, and is beneficial to global geospatial data integration and earth system mode calculation.

(2) Three adjacent rhombic surfaces of the icosahedron are spliced into a combined structure, frequent cross-surface operation among 30 rhombic surfaces is simplified into operation among 10 combined structures, and the processing capacity of large-area geospatial data is improved.

(3) The grid units are identified by adopting quadtree hierarchical codes, the spatial positions of the units are described by adopting a triaxial integer coordinate description unit of a combined structure, and the advantages of the hierarchical codes and the integer coordinate codes are fused, so that a rapid adjacent query method is designed.

Drawings

FIG. 1 is a schematic diagram of a rhombic icosahedron combination structure;

FIG. 2 is an exemplary diagram of a three-axis integer coordinate system of a composite structure;

FIG. 3 is a schematic diagram of hierarchical cutaway and coding;

FIG. 4 is a diagram of the global grid generation result;

FIG. 5 is a diagram illustrating three-axis coordinate conversion into Morton codes;

FIG. 6 is a diagram of a transition seat label for transcoding;

FIG. 7 is a schematic diagram of a unit proximity query;

FIG. 8 is an exemplary diagram of a unit proximity query rule within a composite structure;

FIG. 9 is a diagram showing the transformation of projection coordinates into three-axis integer coordinates;

FIG. 10 is a diagram showing three-axis integer coordinate conversion into a projected coordinate map;

FIG. 11 is a graph showing the comparison of longitude and latitude conversion to coding efficiency;

FIG. 12 is a graph showing the comparison of the transcoding to latitude and longitude efficiency;

FIG. 13 is a comparison of coding proximity query efficiency.

Detailed Description

The invention is further illustrated by the following description of specific embodiments in conjunction with the accompanying drawings:

a method for encoding and operating a rhombic icosahedron global discrete grid system, comprising:

step 1, a rhombic icosahedron combined structure and a unit position description mode are provided: the adjacent upper, middle and lower three diamond surfaces are combined into a structure, as shown in figure 1, the diamond-shaped icosahedron consists of 10 identical combined structures, and the serial numbers are 0 to 9; and establishing a three-axis integer coordinate system IJK in each combined structure, and uniformly describing the spatial positions of all units in the combined structure as shown in fig. 2 without considering the cross-plane problem among three diamond planes in the same combined structure.

And 2, designing a hierarchical segmentation and coding scheme based on the rhombic icosahedron combined structure.

As shown in fig. 3, an initial cell division pattern (fig. 3 (a)) of any combination structure is defined, the upper rhombic-face white hexagonal cell number is 0, the middle rhombic-face gray pentagonal cell number is 1, and the lower rhombic-face blue hexagonal cell number is 2. Adopting different four-hole hierarchical splitting and coding methods (fig. 3 (b)) for different initial units, representing sub-units at different positions by using 0-3 code elements, and adopting a hexagonal quadtree splitting and Morton code coding scheme for the white and blue initial units and the sub-units thereof, so that 3 boundary sub-units are adjacent to the current unit; the gray initial unit and its subunits are hierarchical coded such that 1 subunit is exactly located completely outside the current unit in spatial position. The Morton code Z curve of the layer 2 trellis is shown in FIG. 3 (c), so each cell code is unique and each bit code is associated with a geospatial location. In addition, there is an initial pentagon unit at the north and south vertices of the icosahedron, which is assigned the combination of numbers 2 and 7 and numbered 3, and each layer has only a unique center aligned subunit, and 0 is added to identify the subunits after encoding the units of the layer.

The design level coding rule is as follows, the grid level is n, one Morton code element is added for each layer of the level, the coding form is simple and easy to understand and realize, and the invention is called as a rhombic icosahedron combined structure hexagonal quadtree (Rhombic Triacontahedron Hexagonal Quadtree on Combined Structure, RTHQCS).

Table 1 hierarchical coding scheme

Therefore, global seamless non-stacked mesh subdivision and coding can be realized, as shown in fig. 4, the generation results of the layers 2 and 3 are shown, the unit arrangement in 10 combined structures is identical, and pentagonal units at the south and north vertexes are marked red.

And 3, constructing the interconversion of the hierarchical coding and the triaxial integer coordinates. Specifically comprises the steps 3.1-3.2.

And 3.1, converting the three-axis integer coordinate into code. As shown in FIG. 5, taking the diamond surface on the combined structure as an example, the three axes of the combined structure are first integratedThe numerical coordinates (I, J, K) are subjected to coordinate transformation to obtain diamond transition coordinates (I, J), the transformation method is shown in formula (1), and num=2 ⁿ The number of units on one side of the diamond surface when the grid is layered at present is represented, for example, the num value of the 2 nd grid is 4.

Then adding binary components of the transition coordinates to obtain Morton codes, wherein binary numbers corresponding to the unit transition coordinates i and j are respectively B _i And B _j Morton codes are obtained by adding binary coordinate components as shown in formula (2).

M＝B _i +2B _j (2)

And 3.2, converting the codes into three-axis integer coordinates. As shown in fig. 6, taking a diamond surface on the combined structure as an example, according to the relation between the Morton code and the unit space position, the Morton code is judged layer by layer to perform multiplication and translation of the transition coordinate, and finally the tri-axial integer coordinate can be obtained through inverse transformation, which specifically comprises the steps 3.2.1-3.2.3.

Step 3.2.1, starting from the layer 0 initial grid element, as shown by the black dashed line element in fig. 6, the transition coordinates are i=j=0.

Step 3.2.2, the following operations are performed on each bit Morton, and finally the transition coordinates of the unit can be obtained. First, the transition coordinates i=2i, j=2j of the center subunit of the next layer are calculated. Second, if M [ r ] =0 (r=1, 2, …, n), i.e. the cell is the central subunit, i=i, j=j; if M [ r ] =1, i.e. the cell is the upper left boundary subunit, i=i+1, j=j; if M [ r ] =2, i.e. the cell is the lower boundary subunit, i=i, j=j+1; if M [ r ] =3, i.e. the cell is the lower left boundary subunit, i=i+1, j=j+1.

And 3.2.3, converting the transition coordinates into three-axis integer coordinates, wherein the process is the inverse process of the formula (1).

It should be noted that, the processing manner of the middle diamond surface and the lower diamond surface is similar to that of the upper diamond surface, and will not be described herein.

And 4, designing a coding operation method for the RTHQCS. Specifically comprises the steps 4.1-4.2.

And 4.1, designing a hierarchical coding operation method. The method comprises the steps of searching parent unit codes and searching child unit codes, wherein the hierarchy subdivision rule of the method is obtained, and only the tail bits of the codes are needed to be removed when the parent units are searched; when searching the subunit, the symbol set {0,1,2,3} is added to the last bit of the code according to the space position of the subunit and the unit.

And 4.2, designing a neighboring coding operation method. As shown in fig. 7, the neighbor operation of the hexagonal cell is divided into three cases: the hexagons and their neighboring cells are all located inside the belonging composite structure (referred to herein as internal cells), the hexagons are located on the boundary of the belonging composite structure (referred to herein as boundary cells), the hexagons are located inside the belonging composite structure but some of their neighboring cells are affiliated to the neighboring composite structure (referred to herein as half-boundary cells). The three units of proximity operation specifically includes steps 4.2.1-4.2.3.

Step 4.2.1, designing a neighbor operation method of the internal unit. As shown in the gray cell in fig. 7, the proximity operation is performed by simply adding or subtracting 1 to or from the three-axis integer coordinate without considering the cross-plane problem, and is divided into the following three cases: when the unit is located between the I axis and the K axis, the coordinates are (I, 0, K), and the rule of the adjacent operation is as shown in fig. 8 (a); when the unit is located on the I axis, the coordinates are (I, 0), and the rule of the neighbor operation is as shown in fig. 8 (b); when the unit is located between the I axis and the J axis, the coordinates are (I, J, 0), and the rule of the neighbor operation is as shown in fig. 8 (c).

Step 4.2.2, designing a neighboring operation method of the boundary unit. As shown in the green cell in fig. 7, adjacent cells in the belonging combined structure are first searched by the method of fig. 8, and other adjacent cells are located on the adjacent combined structure, and a cross-plane adjacent search method needs to be studied. The combined structures with the numbers of 0 to 4 are classified into a first class, the first class is marked as S1, the combined structures with the numbers of 5 to 9 are correspondingly marked as S2, and the upper diamond surface, the middle diamond surface and the lower diamond surface of the combined structures are respectively represented by (0), (1) and (2). The invention discovers that the coding cross-plane operation of the units between the adjacent combined structures has a certain rule, and particularly refers to a code element replacement rule between the coding of the boundary unit and the coding of the upper half boundary unit of the other combined structure. The unit positions can be divided into 3 types of Morton symbol substitution methods, and Morton symbols before conversion are represented by M r (r=1, 2, …, n), and Morton symbols after conversion are represented by M' r, as shown in Table 2.

Table 2 symbol substitution method for cross-combining structure proximity query

And 4.2.3, designing a proximity operation method of the half boundary unit. As shown in the blue unit in fig. 7, the adjacent units in the associated combination structure are searched first by the method of fig. 8, and then other adjacent units on the adjacent combination structure are found by the inverse process of the symbol substitution, so that the inverse transformation is easy to be implemented, and the invention is not repeated.

And 5, constructing the mutual conversion of the geographic coordinates and the grid codes, and ensuring that the common geographic space data are smoothly integrated into the grid system. Specifically comprises a step 5.1 and a step 5.2.

And 5.1, converting the geographic coordinates into codes. Firstly, projecting a geographic coordinate to an icosahedron diamond-shaped surface to obtain a projection coordinate (x, y); then converting the projection coordinates into three-axis integer coordinates (I, J, K) according to the grid subdivision geometric characteristics; finally, the three-axis integer coordinates are converted into codes by the method of step 3.1. Taking a diamond surface on a combined structure as an example, the specific process of converting projection coordinates into triaxial integer coordinates is as follows:

firstly, knowing the projection coordinates (x, y) of any point P, as shown in fig. 9, under the projection coordinate system, a straight line equation of a right lower diamond side, namely a K axis, is y=tan (58.285 °) x-cos (31.715 °) L, and a straight line equation of an I axis is y= -cos (31.715 °) L, wherein L is a diamond side length; respectively calculating Euclidean distances d between points P, K axis and I axis _I And d _K Then the point P is decomposed into floating point coordinates l on the coordinate axes I and K by utilizing the geometric relation _I And l _K The calculation formula is as follows, whereinFor the spacing of adjacent hexagonal cell centers in the directions of the K-axis and the J-axis +.>Is the distance between adjacent hexagonal cell centers along the I-axis direction, and the two satisfy the relation:

second, m _I And m _K Respectively is l _I And l _K As a result of rounding, the fractional part of point P on the I-axis is r _I ＝l _I -m _I The fractional part on the K-axis is r _K ＝l _K -m _K R corresponding to three points A, B and C in FIG. 9 _K The values are respectivelyThe integer coordinates of the unit can thus be calculated by the formula (4), wherein +_ is obtainable according to the triangle sine theorem>

Finally, as shown by the red dashed line elements in fig. 9, these boundary elements are attached to adjacent composite structures and require cross-plane processing. Firstly, subtracting 1 from the I coordinate, namely translating to a half boundary unit adjacent to the I coordinate, then converting the three-axis integer coordinate of the adjacent unit into codes, and finally obtaining the codes of the boundary units through the inverse process of code element replacement shown in the table 2.

And 5.2, converting the codes into geographic coordinates. Firstly, converting codes into three-axis integer coordinates (I, J, K) by a method of the step 3.2; then calculating diamond plane projection coordinates (x, y); finally, the geographical coordinates (B, L) of the cell centers are obtained through back projection calculation. Taking a diamond surface on a combined structure as an example, the specific process of converting the three-axis integer coordinate into the unit center projection coordinate is as follows:

as shown in FIG. 10, if the coordinates corresponding to the arbitrary point P are (I, 0, K), the coordinates of the upper diamond-shaped surface center unit O areAnd the grid center is the origin of the projection coordinate system. First, the coordinate difference between the point P and the center unit O is calculated by equation (5):

the projected coordinates of the point P are then calculated by the geometric relationship as shown in equation (6).

It should be noted that, the processing manner of the middle diamond surface and the lower diamond surface is similar to that of the upper diamond surface, and will not be described herein.

And 6, performing quality evaluation and comparison on the rhombic icosahedron global discrete grid system established by the invention. The prior rhombic icosahedron grid system has no hierarchical coding scheme, therefore, the invention selects HLQT and HHOT with higher efficiency in the regular icosahedron four-hole hexagonal grid as comparison objects, traverses the world by 0.2 degree longitude and latitude intervals, selects 162 ten thousand points in total, selects 13-20 high-level grids, respectively tests the average time of three operations of converting longitude and latitude coordinates into coding, converting the coding into longitude and latitude coordinates and coding adjacent query under similar resolution. Each operation was run 10 times and averaged in microseconds (mus), with the efficiency ratio being the ratio of HLQT or HHUT for use with respect to RTHQBS for use. All programs were compiled to Release version and tested in a compatible machine (hardware configuration: interCore i5-10400F [email protected],16G RAM,KIOXIA-EXCERIA480G SSD; operating system: windows 10x64 Enterprise LTSC; development tool: microsoft Visual C ++ Enterprise 2022).

Because the projection methods adopted by different DGGS schemes are different, the influence of different projection modes on the conversion efficiency is different in the process of interconversion between longitude and latitude and codes. Therefore, in order to more fairly compare various coding schemes, the time of forward projection and backward projection is eliminated, and only the interconversion time of the projection coordinates of the polyhedral surface and the codes is considered.

As shown in fig. 11, the efficiency of the transformation of the longitude and latitude coordinates of the RTHQCS into codes is highest at each level (fig. 11 (a)), and the efficiency advantage gradually increases with the level increase (fig. 11 (b)), which can reach 27.01 times of HLQT and 4.29 times of HHUT. This benefits from the fast conversion algorithm of the three-axis integer coordinates to the quaternary Morton codes of the present invention, and the conversion algorithm is hardly affected by the grid level. The HLQT consists of 12 vertex tiles and 20 surface tiles, the structure is more complex, in addition, the HLQT and HHUT all need to perform coding addition operation one symbol by one symbol to obtain final coding, and the coding length increases along with the increase of the layers, so that the calculated amount of the coding addition increases.

As shown in fig. 12, the efficiency of transcoding the RTHQCS scheme to latitude and longitude is highest at each level (fig. 12 (a)), and can reach 46.86 times HLQT and 2.36 times HHUT as the level increases with the efficiency advantage and increases gradually (fig. 12 (b)). The method is characterized in that the code element of the scheme only has 0,1,2 and 3, the structure is simple, the judgment condition is few, the code element does not need to participate in calculation when being 0, and the projection plane coordinate can be obtained through simple geometric relation calculation after the code is converted into the triaxial integer coordinate; each layer of code element of HLQT needs to participate in calculation, and the first code element is complex, thus consuming a great deal of processing time; although the HHOT does not need to participate in calculation when the code element is 0, only the Cartesian coordinates of the tile can be obtained after scaling and accumulating by the complex number codes, and multiple judgment is needed to calculate the projection plane coordinates.

As shown in fig. 13, the encoding proximity query efficiency of the RTHQCS scheme is highest at each level (fig. 13 (a)), and the efficiency ratio gradually increases as the level increases (fig. 13 (b)), which can reach 40.97 times HLQT, 4.35 times HHUT. The number of the faces of the rhombic thirty faces is larger than that of the regular icosahedron, accordingly, more cross-face operation is involved in adjacent query, so that the adjacent query efficiency is reduced, but the adjacent three rhombic faces are combined into a combined structure and a three-axis integer coordinate system is established, the problem of adjacent query of the cross-rhombic faces in the combined structure is solved, the cross-face problem between the adjacent combined structures can be directly realized through code element replacement, the complexity of the adjacent query is mainly related to the efficiency of code conversion into the three-axis integer coordinate, and therefore, the adjacent query efficiency of the scheme is reduced the slowest with the increase of layers. The adjacent query of the HLQT is realized through the coding addition operation with 6 direction vectors, and the efficiency is lowest; the HHUT scheme part of the neighboring cells can be obtained directly by changing the last bit coding, but other cells also have to resort to coding addition. Both of these two comparison schemes require searching the code addition lookup table symbol by symbol, so as the trellis level increases, the code length increases, and the computational efficiency becomes slower.

On the basis of the embodiment, the invention also provides a coding and operation device of the rhombic icosahedron global discrete grid system, which comprises:

the first construction module is used for dividing the rhombic icosahedron into 10 combined structures, the serial numbers of the combined structures are 0-9 in sequence, and each combined structure is formed by combining an upper diamond surface, a middle diamond surface and a lower diamond surface which are adjacent; establishing a triaxial integer coordinate system in each combined structure;

the hierarchical splitting and encoding module is used for constructing an initial unit for each combined structure, wherein the initial unit comprises an upper rhombic-surface hexagonal unit, a middle rhombic-surface pentagonal unit and a lower rhombic-surface hexagonal unit, the upper rhombic-surface hexagonal unit is numbered to be 0, the middle rhombic-surface pentagonal unit is numbered to be 1, and the lower rhombic-surface hexagonal unit is numbered to be 2; in addition, the initial pentagonal units at the south and north vertexes of the rhombic icosahedron are numbered 3 and respectively assigned to a 2-7 combined structure; different four-hole hierarchical splitting and coding methods are adopted for different initial units;

the second construction module is used for constructing the interconversion between the hierarchical coding and the triaxial integer coordinates;

the coding operation design module is used for designing a coding operation method;

and the third construction module is used for constructing the mutual conversion between the geographic coordinates and the grid codes so as to integrate the geographic space data into the grid system and obtain the rhombic icosahedron global discrete grid system.

Further, the four-hole hierarchical subdivision and coding method for different initial units comprises the following steps: for the initial pentagon unit of the north and south vertexes, only a unique sub-unit with aligned centers is identified by a code element 0 after each layering.

Further, the four-hole hierarchical subdivision and coding method for different initial units comprises the following steps:

performing four-hole hierarchical segmentation and coding on other initial units except for the initial pentagon unit of the north-south vertex, representing subunits at different positions by 0-3 code elements, wherein 0 is used for identifying a central subunit, 1-2 is used for identifying a boundary subunit, 3 is used for identifying a boundary subunit or an external subunit, and the hexagonal units of the upper diamond surface and the hexagonal units of the lower diamond surface and the two subunits adopt a hexagonal quadtree segmentation and Morton code coding scheme, so that 3 boundary subunits are adjacent to the current unit; and carrying out hierarchical subdivision coding on the pentagon unit and the sub-units of the pentagon unit with the middle diamond surface, so that 1 sub-unit is completely positioned outside the current unit in the space position.

Further, the second construction module is specifically configured to perform the following steps:

step 3.1, converting the triaxial integer coordinates into codes, and for the diamond-shaped surface on the combined structure, firstly, carrying out coordinate transformation on the triaxial integer coordinates (I, J, K) of the combined structure according to the formula (1) to obtain diamond-shaped surface transition coordinates (I, J):

wherein num=2 ⁿ Representing the number of units on one side of the diamond-shaped surface when the current grid is layered;

the binary components of the transition coordinates are then added according to equation (2) to obtain Morton codes:

M＝B _i +2B _j (2)

wherein B is _i And B _j Binary numbers corresponding to the transition coordinates i and j of the units are respectively represented;

and 3.2, converting the codes into triaxial integer coordinates, judging Morton codes layer by layer according to the relation between the Morton codes and the unit space positions for the diamond surface on the combined structure, multiplying and translating the transition coordinates, and finally obtaining the triaxial integer coordinates through inverse transformation.

Further, the step 3.2 includes:

step 3.2.1, setting the transition coordinates of the initial grid unit of the 0 th layer to be (i, j) = (0, 0);

step 3.2.2, executing the following operations on each Morton code to finally obtain the transition coordinates of the unit: calculating transition coordinates i=2i, j=2j of a central subunit of the next layer, and if the r-th Morton symbol Mr of the subunit is 0, i.e. the unit is the central subunit, i=i, j=j; if M [ r ] =1, i.e. the cell is the upper left boundary subunit, i=i+1, j=j; if M [ r ] =2, i.e. the cell is the lower boundary subunit, i=i, j=j+1; if M [ r ] =3, i.e. the cell is the lower left boundary subunit, i=i+1, j=j+1;

and 3.2.3, converting the transition coordinates into three-axis integer coordinates, wherein the process is the inverse process of the formula (1).

Further, the encoding operation design module is specifically configured to execute the following steps:

step 4.1, designing a hierarchical coding operation method, namely searching for a parent unit code and searching for a child unit code, and removing the tail bit of the code when searching for the parent unit code; when searching for the subunit codes, adding a code element set {0,1,2,3} in the tail bits of the codes according to the space positions of the subunit and the father unit;

step 4.2, designing a neighboring encoding operation method, wherein the neighboring operation of the hexagonal unit is divided into three cases: the hexagons and adjacent units thereof are all positioned inside the affiliated combined structure, the hexagons are positioned on the boundary of the affiliated combined structure, the hexagons are positioned inside the affiliated combined structure, but part of adjacent units thereof are affiliated to the adjacent combined structure, the units corresponding to the first case are called as internal units, the units corresponding to the second case are called as boundary units, and the units corresponding to the third case are called as half-boundary units.

Further, the step 4.2 includes:

step 4.2.1, completing the adjacent operation of the internal unit by adding or subtracting 1 to the three-axis integer coordinate;

step 4.2.2, for boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, wherein other adjacent units are positioned on the adjacent combined structure, classifying the combined structure with the number of 0-4 into a first class, marking the combined structure as S1, correspondingly marking the combined structure with the number of 5-9 as S2, respectively using (0), (1) and (2) to represent upper, middle and lower diamond surfaces of the combined structure, and inquiring the adjacent combined structure half boundary units by the boundary units according to different unit positions to obtain 3 Morton code element replacing methods as follows:

wherein M r represents Morton symbols of the unit before transformation, and M' r represents Morton symbols after transformation;

step 4.2.3, for half boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, and then finding other adjacent units on the adjacent combined structure by the reverse process of symbol substitution.

Further, the third construction module is specifically configured to perform the following steps:

step 5.1, converting the geographic coordinates into codes: firstly, projecting a geographic coordinate to an icosahedron diamond-shaped surface to obtain a projection coordinate (x, y); then converting the projection coordinates into three-axis integer coordinates (I, J, K) according to the grid subdivision geometric characteristics; finally, converting the three-axis integer coordinates into codes by the method of the step 3.1;

step 5.2, converting the codes to geographic coordinates: firstly, converting codes into three-axis integer coordinates (I, J, K) by a method of the step 3.2; then calculating diamond plane projection coordinates (x, y); finally, the geographical coordinates (B, L) of the cell centers are obtained through back projection calculation.

Further, the method further comprises the following steps:

the quality evaluation and comparison module is used for evaluating and comparing the quality of the established rhombic icosahedron global discrete grid system.

In summary, the invention proposes to splice three adjacent rhombic surfaces of the icosahedron into a combined structure, and simplify the frequent cross-surface operation among 30 rhombic surfaces into the operation among 10 combined structures, thereby improving the processing capacity of large-area geospatial data. The grid units are identified by adopting quadtree hierarchical codes, the spatial positions of the units are described by adopting a triaxial integer coordinate description unit of a combined structure, and the advantages of the hierarchical codes and the integer coordinate codes are fused, so that a rapid adjacent query method is designed.

The foregoing is merely illustrative of the preferred embodiments of this invention, and it will be appreciated by those skilled in the art that changes and modifications may be made without departing from the principles of this invention, and it is intended to cover such modifications and changes as fall within the true scope of the invention.

Claims (10)

1. The coding and operation method of the rhombic icosahedron global discrete grid system is characterized by comprising the following steps of:

step 1, dividing a rhombic icosahedron into 10 combined structures, wherein the serial numbers of the combined structures are 0-9, and each combined structure is formed by combining adjacent upper, middle and lower rhombic faces; establishing a triaxial integer coordinate system in each combined structure;

step 2, for each combined structure, constructing an initial unit comprising an upper diamond-shaped hexagonal unit, a middle diamond-shaped pentagonal unit and a lower diamond-shaped hexagonal unit, wherein the number of the upper diamond-shaped hexagonal unit is 0, the number of the middle diamond-shaped pentagonal unit is 1, and the number of the lower diamond-shaped hexagonal unit is 2; in addition, the initial pentagonal units at the south and north vertexes of the rhombic icosahedron are numbered 3 and respectively assigned to a 2-7 combined structure; different four-hole hierarchical splitting and coding methods are adopted for different initial units;

step 3, constructing the interconversion of the hierarchical coding and the triaxial integer coordinates;

step 4, designing a coding operation method;

and 5, constructing the mutual conversion of the geographic coordinates and the grid codes so as to integrate the geographic space data into the grid system and obtain the rhombic icosahedron global discrete grid system.

2. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 1, wherein the method for adopting different four-hole hierarchical subdivision and encoding for different initial units comprises: for the initial pentagon unit of the north and south vertexes, only a unique sub-unit with aligned centers is identified by a code element 0 after each layering.

3. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 2, wherein the method for adopting different four-hole hierarchical subdivision and encoding for different initial units comprises:

performing four-hole hierarchical segmentation and coding on other initial units except for the initial pentagon unit of the north-south vertex, representing subunits at different positions by 0-3 code elements, wherein 0 is used for identifying a central subunit, 1-2 is used for identifying a boundary subunit, 3 is used for identifying a boundary subunit or an external subunit, and the hexagonal units of the upper diamond surface and the hexagonal units of the lower diamond surface and the two subunits adopt a hexagonal quadtree segmentation and Morton code coding scheme, so that 3 boundary subunits are adjacent to the current unit; and carrying out hierarchical subdivision coding on the pentagon unit and the sub-units of the pentagon unit with the middle diamond surface, so that 1 sub-unit is completely positioned outside the current unit in the space position.

4. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 1, wherein said step 3 comprises:

step 3.1, converting the triaxial integer coordinates into codes, and for the diamond-shaped surface on the combined structure, firstly, carrying out coordinate transformation on the triaxial integer coordinates (I, J, K) of the combined structure according to the formula (1) to obtain diamond-shaped surface transition coordinates (I, J):

wherein num=2 ⁿ Representing the number of units on one side of the diamond-shaped surface when the current grid is layered;

the binary components of the transition coordinates are then added according to equation (2) to obtain Morton codes:

M＝B _i +2B _j (2)

wherein B is _i And B _j Binary numbers corresponding to the transition coordinates i and j of the units are respectively represented;

and 3.2, converting the codes into triaxial integer coordinates, judging Morton codes layer by layer according to the relation between the Morton codes and the unit space positions for the diamond surface on the combined structure, multiplying and translating the transition coordinates, and finally obtaining the triaxial integer coordinates through inverse transformation.

5. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 4, wherein said step 3.2 comprises:

step 3.2.1, setting the transition coordinates of the initial grid unit of the 0 th layer to be (i, j) = (0, 0);

step 3.2.2, executing the following operations on each Morton code to finally obtain the transition coordinates of the unit: calculating transition coordinates i=2i, j=2j of a central subunit of the next layer, and if the r-th Morton symbol Mr of the subunit is 0, i.e. the unit is the central subunit, i=i, j=j; if M [ r ] =1, i.e. the cell is the upper left boundary subunit, i=i+1, j=j; if M [ r ] =2, i.e. the cell is the lower boundary subunit, i=i, j=j+1; if M [ r ] =3, i.e. the cell is the lower left boundary subunit, i=i+1, j=j+1;

and 3.2.3, converting the transition coordinates into three-axis integer coordinates, wherein the process is the inverse process of the formula (1).

6. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 1, wherein said step 4 comprises:

step 4.1, designing a hierarchical coding operation method, namely searching for a parent unit code and searching for a child unit code, and removing the tail bit of the code when searching for the parent unit code; when searching for the subunit codes, adding a code element set {0,1,2,3} in the tail bits of the codes according to the space positions of the subunit and the father unit;

step 4.2, designing a neighboring encoding operation method, wherein the neighboring operation of the hexagonal unit is divided into three cases: the hexagons and adjacent units thereof are all positioned inside the affiliated combined structure, the hexagons are positioned on the boundary of the affiliated combined structure, the hexagons are positioned inside the affiliated combined structure, but part of adjacent units thereof are affiliated to the adjacent combined structure, the units corresponding to the first case are called as internal units, the units corresponding to the second case are called as boundary units, and the units corresponding to the third case are called as half-boundary units.

7. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 6, wherein said step 4.2 comprises:

step 4.2.1, completing the adjacent operation of the internal unit by adding or subtracting 1 to the three-axis integer coordinate;

step 4.2.2, for boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, wherein other adjacent units are positioned on the adjacent combined structure, classifying the combined structure with the number of 0-4 into a first class, marking the combined structure as S1, correspondingly marking the combined structure with the number of 5-9 as S2, respectively using (0), (1) and (2) to represent upper, middle and lower diamond surfaces of the combined structure, and inquiring the adjacent combined structure half boundary units by the boundary units according to different unit positions to obtain 3 Morton code element replacing methods as follows:

wherein M r represents Morton symbols of the unit before transformation, and M' r represents Morton symbols after transformation;

step 4.2.3, for half boundary units, searching adjacent units in the combined structure by the method of step 4.2.1, and then finding other adjacent units on the adjacent combined structure by the reverse process of symbol substitution.

8. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 4, wherein said step 5 comprises:

step 5.1, converting the geographic coordinates into codes: firstly, projecting a geographic coordinate to an icosahedron diamond-shaped surface to obtain a projection coordinate (x, y); then converting the projection coordinates into three-axis integer coordinates (I, J, K) according to the grid subdivision geometric characteristics; finally, converting the three-axis integer coordinates into codes by the method of the step 3.1;

step 5.2, converting the codes to geographic coordinates: firstly, converting codes into three-axis integer coordinates (I, J, K) by a method of the step 3.2; then calculating diamond plane projection coordinates (x, y); finally, the geographical coordinates (B, L) of the cell centers are obtained through back projection calculation.

9. The method for encoding and computing a rhombic icosahedron global discrete grid system according to claim 1, further comprising, after said step 5:

and 6, performing quality evaluation and comparison on the established rhombic icosahedron global discrete grid system.

10. A coding and computing device of a rhombic icosahedron global discrete grid system, comprising:

the first construction module is used for dividing the rhombic icosahedron into 10 combined structures, the serial numbers of the combined structures are 0-9 in sequence, and each combined structure is formed by combining an upper diamond surface, a middle diamond surface and a lower diamond surface which are adjacent; establishing a triaxial integer coordinate system in each combined structure;

the hierarchical splitting and encoding module is used for constructing an initial unit for each combined structure, wherein the initial unit comprises an upper rhombic-surface hexagonal unit, a middle rhombic-surface pentagonal unit and a lower rhombic-surface hexagonal unit, the upper rhombic-surface hexagonal unit is numbered to be 0, the middle rhombic-surface pentagonal unit is numbered to be 1, and the lower rhombic-surface hexagonal unit is numbered to be 2; in addition, the initial pentagonal units at the south and north vertexes of the rhombic icosahedron are numbered 3 and respectively assigned to a 2-7 combined structure; different four-hole hierarchical splitting and coding methods are adopted for different initial units;

the second construction module is used for constructing the interconversion between the hierarchical coding and the triaxial integer coordinates;

the coding operation design module is used for designing a coding operation method;

and the third construction module is used for constructing the mutual conversion between the geographic coordinates and the grid codes so as to integrate the geographic space data into the grid system and obtain the rhombic icosahedron global discrete grid system.