Precise hexagonal pixel modeling and an easy-sharing storage scheme for remote sensing images based on discrete global grid system

ABSTRACT

Discrete Global Grid Systems (DGGSs) are an emerging Earth reference model that support the integration and analysis of remote sensing (RS) data. Grid modeling, sampling, quantization, and storage are the key points and difficulties of DGGS. The Icosahedral Snyder Equal-area Aperture-4 Hexagonal DGGS was introduced as the basic framework to improve the geospace sampling efficiency. We presented an approach for precise hexagonal pixel modeling and an easy-sharing storage scheme compatible with open-standard formats for RS images based on the DGGS. Firstly, the proposed interpolation is computed by overlaps between quadrilateral and hexagonal pixels. The hexagonal grids were then mapped onto the icosahedral surface, and a strict correspondence between hexagonal and rectangular pixels is established. The open-standard formats were used to accurately store the hexagonal attribute values and metadata. Finally, a multiscale hexagonal aggregation algorithm based on the dataset was designed. Experiments showed that the proposed modeling methodology had a higher accuracy and smaller errors. The hexagonal data can be stored as regular rectangles with a fixed pattern in any standard format. This storage scheme was more conducive to data processing and sharing compared to SMOS, which was expected to promote hexagonal DGGS in RS data organization, processing and sharing.

KEYWORDS:

• Discrete global grid systems
• hexagonal pixels
• quadrilateral pixels
• mathematical model
• remote sensing data
• GeoTIFF

1. Introduction

The integration and analysis of large amounts of multisource, multiscale, and multidimensional remote sensing (RS) data have become a common application requirement because of the rapid development of Earth observation technologies and continuous growth in available geospatial data (Guo et al. 2017). There are data cubes and spatio-temporal databases (e.g. GeoMesah or MoblityDB) that support storage and analysis in both the temporal and spatial dimensions of data, such as the Australian Geoscience Data Cube (AGDC) (Lewis et al. 2017). Conventional two-dimensional (2-D) planar data models based on local map projections hardly support global-scale analyses (Mahdavi-Amiri, Alderson, and Samavati 2015). There is also area deterioration at the poles in three-dimensional (3-D) latitude and longitude grids. Discrete Global Grid Systems (DGGSs) have been proposed as the fundamental framework for managing geospatial data (Sahr, White, and Kimerling 2003) and a potential data model for supporting computational infrastructure (Robertson et al. 2020). It recursively discretizes Earth’s surface using a specialized method to form a seamless and non-overlapping hierarchy (Open Geospatial Consortium 2017). More critically, DGGSs are better suited for handling large-scale issues, and are expected to offer solutions for new digital Earth data processing and analysis (Mahdavi-Amiri, Alderson, and Samavati 2015).

There are three regular subdivision grids: triangles, quadrilaterals, and hexagons (Kiselman 2022; Klette and Rosenfeld 2004), and some semi-regular grids (Nagy 2022). In the applications of regular grids, hexagons have more advantages for geospatial data and image processing because of their good geometric properties compared with triangles and quadrilaterals (Middleton and Sivaswamy 2005), such as consistent connectivity (Liao et al. 2020), higher sampling efficiency (Mersereau 1979), and higher angular resolution (Wang et al. 2020).

The information present in the RS data must be accurately converted to a hexagonal grid representation to fully exploit the benefits of hexagonal DGGS (Robertson et al. 2020). This is an essential requirement for the subsequent processing of grid data, cell aggregation, and spatial query access. Although the Open Geospatial Consortium (OGC) does not impose any restrictions on the techniques used to convert data to DGGS, an accepted quantization strategy is crucial for accuracy control (Li and Stefanakis 2020). Ma et al. (2021) demonstrated that bilinear interpolation preserves most structural features and information. However, bilinear interpolation may result in an invalid value because there are insufficient appropriate interpolation inputs close to the boundary of the image block, which affects the completeness of the global data (Li, McGrath, and Stefanakis 2021).

According to several significant studies on hexagonal DGGS, wavelet transformation with other technologies has been introduced (Mahdavi-Amiri, Alderson, and Samavati 2016) to manage multiscale geospatial data and achieve effective network transmission using an aperture-3 hexagonal grid system. An online geospatial data analysis platform (Global Grid Systems 2023) was also created based on an aperture-3 hexagonal hierarchy, providing customers with features, such as integrated management and real-time processing of multisource data. Uber integrated a large amount of taxi data using an aperture-7 hexagonal grid system to achieve reasonable scheduling and management of vehicles (Uber 2023).

Although the aperture-3 hexagonal hierarchy can achieve smooth transitions between levels, the orientation of the grids rotates by 30 degree in successive resolution, complicating the coding scheme. The aperture-7 subdivision offers effective hierarchical operations for successive levels; however, it is difficult to select an appropriate cell when aggregating data (Zhou et al. 2022). Aperture-4 hexagonal DGGS has a fixed orientation, facilitating the retrieval and aggregation of cells (Guo, Goodchild, and Annoni 2020). This strategy can easily be combined with extensive operations. The European Space Agency (ESA) collects soil moisture and ocean salinity (SMOS) data on aperture-4 hexagonal grids and makes the data available to users in the form of an approximate hexagon comprising several dozen square pixels (Martin Suess 2004). Despite having great efficiency in visualization, this not only increases the volume of data, but also causes severe distortion as the latitude increases.

Some researchers have used ASCII (Sousa and Leitão 2018), GeoJSON (Li, McGrath, and Stefanakis 2021), or other file formats to store hexagonal data. This scheme is particularly convenient for time-series streaming, but appears unsuitable for hexagonal grids with specific geographic positions and ranges. The topological structure of the cells must be explicitly stored, thus increasing the value associated with a single hexagon from 1 to 13 (adding the coordinates of each vertex). Currently, there is no standard file format for hexagonal data; as hexagons cannot be arranged in the same manner as squares, existing storage methods cannot be utilized directly. The private file formats often used by scholars are not straightforward for direct interpretation using commercial software, making sharing with others challenging.

Based on this analysis, we established a mathematical model to represent hexagonal RS data and introduced a storage scheme using the standard GeoTIFF file for encoding, and processing geographic data. It is not only convenient for analyzing data in general software, but also takes full advantage of the hexagon. The remainder of this paper is organized as follows. Section 2 introduces the basic idea of this study. Section 3 provides the methodology of precise hexagonal pixel modeling and an organization scheme for the hexagonal RS data. Section 4 describes the experimental process and presents the results. Finally, Section 5 concludes the study and provides suggestions for future research.

2. Basic ideas

This section discusses the mathematical model for representing hexagonal pixel information from RS data using Icosahedral Snyder Equal-area Aperture-4 Hexagonal (ISEA4H) DGGS as well as the key process of organizing hexagonal pixel data using an open standard format. The basic concept is illustrated in Figure 1.

Figure 1. Basic concepts of this study.

Traditionally, RS datasets are mainly referenced in a continuous space via the projected Cartesian space or represented in latitude and longitude. However, the resolution of datasets often varies with different surface regions; latitude and longitude grids suffer varying degrees of deformation from the equator to the poles (Kelly and Šavrič 2021), leading to redundancy in data representation (Figure 1(a)). Hence, Cartesian grid-based data representations are considered antiquated, imprecise, and unsuitable for global data analysis (Baumann 2021). In particular, for RS data modeling and organization, grids with consistency and equal areas are imperative requirements to ensure the accuracy of information collection. Therefore, this study employed hexagonal grid systems as a basis for organizing data to realize uniform sampling of Earth's space and utilized a spherical quantitating method that considers the quadrilateral pixel weight (Figure 1(b)).

Despite studies have shown that the hexagon is the optimal cell for data sampling. By contrast, most graphics hardware do not currently support rendering of hexagonal pixel images, mostly images on square and rectangular grids are supported. They record geographic phenomena by establishing a direct connection between the location information and descriptive property values. Additionally, given the structural characteristics of a rectangle or square, they can be arranged in a reasonably simple and regular manner. Topological connection relations between cells or vertices can be stored implicitly. The value at a certain position can be obtained using a simple and efficient calculation method. GeoTIFF, which is frequently used in RS and GIS, is the standard format for storing square raster data: it is characterized by scalability, interoperability, and ease of shareability.

It is difficult to establish a direct connection between the hexagon and rectangle on the sphere, as well as break the barrier between the hexagonal and rectangular pixels in the RS data. Therefore, the proposed concept can be described as follows. First, we used an icosahedron Snyder equal-area projection to achieve mapping between spherical and polyhedral hexagons (Figure 1(c)). We then considered the unfolding of the icosahedron, flattening the grid onto a 2-D plane that could be easily manipulated and stored (Figure 1 (d)). Next, equal-area and unique correspondences were constructed between the regular hexagon and rectangle on the polyhedral surface. The scheme allowed for a 10-region division of Earth’s surface into rhombic units. Hexagonal pixels were encoded into a GeoTIFF container in the form of rectangular pixels, where each rectangle represents a hexagon, along with the additional parameters and metadata required to reconstruct the hexagonal RS data (Figure 1(e)).

Furthermore, this study created a hexagonal aggregation scheme to manage global multiscale hexagonal RS dataset more efficiently. Owing to the unique link between hexagons and rectangles, rectangular operations were treated as equivalent to hexagonal operations, such that it was easy to quickly obtain spherical hexagonal pixels at various resolutions using the parameters encapsulated in GeoTIFF (Figure 1(f)).

3. Methodology

3.1. Precise mathematical modeling for hexagonal pixel

3.1.1. Grid level determination

The appropriate grid level should first be determined based on the RS spatial resolution for modeling and sampling. The pixel sizes in the RS images and the hexagonal grid with different levels of aperture-4 were both fixed (Thompson et al. 2022), making it difficult to ensure a perfect match between them. The grid level can be determined using equations (1) to be as close to the original image and minimize information loss. (1) $\mathrm{Min}(S_q-S_h)\mathrm{AND}(S_h<S_q)$(1) where$S_h$is the area of a hexagon at level$n$and$S_q$is the area of the RS image pixel.

Taking the commonly used RS datasets as an example, the corresponding grid levels in the different types of images were calculated, as shown in Table 1.

Table 1. Grid levels corresponding to commonly used high-, medium-, and low-resolution images.

Data Type	Remote sensing Data	Spatial resolution (m)	Corresponding grid level
Low spatial resolution	MODIS	250–1000	13
	NOAA/AVHRR	1090	13
Medium spatial resolution	Landsat8	30	18
	Sentinel-2	10	19
High spatial resolution	GF-2	1–4	22
	World View-2	0.46-1	22

3.1.2. Hexagonal pixel modeling

There is currently a wealth of exploration in the theory (Mersereau 1979; Nagy 2015), applications (Brimkov and Barneva 2005; Deutsch 1972; Lukić and Nagy 2019; Matej, Herman, and Vardi 1998), and precision evaluation (Ma et al. 2021) of image processing on the planar hexagonal grid. Related studies have also already been conducted on implementing spherical hexagonal grid to represent RS or terrain data (Sun, Zhao, and Chen 2007; Teanby 2006). The main idea is to project a spherical cell onto the projection plane where the image is located, and then assign attribute values to the cell using linear interpolation. However, there are various projection coordinate systems for multisource images that not only result in time-consuming calculations, but produce varying degrees of distortion, which ultimately affect the accuracy.

Considering that the essence of DGGS is to tessellate the Earth’s surface instead of the plane, it is advisable to complete the modeling on the sphere as much as possible. In addition, a single hexagon may be associated with multiple quadrilateral pixels because they are spatially different types of pixels (Figure 2). Each quadrilateral pixel had a varying degree of impact on the hexagonal pixel value. Consequently, the spatial relationship (e.g. intersection or separation) between the two types of pixels must be properly considered.

Figure 2. Intersections of hexagonal and quadrilateral pixels at different latitudes.

The majority of datasets still employ a traditional model based on map projection owing to the limitations of operating habits and product compatibility. The spatial datum must be aligned under a unified coordinate system before modeling. The spatial intersection relationship can be established by aligning the planar RS image with the grid system (Figure 3). We assumed that the sequence of the hexagonal and quadrilateral vertices was $h_1{h}_2{h}_3{h}_4{h}_5{h}_6,q_1{q}_2{q}_3{q}_4$.

Figure 3. Spherical sampling schematic diagram. (a) Intersection of a quadrilateral and a hexagon. (b) Hexagon is contained in a quadrilateral.

The hexagonal value is equivalent to the quadrilateral value when it is contained within a quadrilateral pixel because every part of the entire hexagon belongs to the same quadrilateral pixel. The hexagonal value is calculated based on the weighted sum of several quadrilateral pixel values when they intersect. Weights were achieved by the overlap area between the target hexagonal pixel and quadrilateral pixels involved in the conversion process. As the intersection area is larger, the weight is higher, which means that the hexagonal value is more affected by the quadrilateral. Suppose that a hexagon intersects with$M$quadrilateral pixels, the hexagonal value,$A$, is calculated as follows: (2) $A={\displaystyle \frac{{\sum }_{i=1}^M(q_i\times \mathrm{Valu}{e}_i)}{{\sum }_{i=1}^M{q}_i}}$(2) where $q_i$is the ratio of the area of the intersecting cell to the quadrilateral pixel and $\mathrm{Valu}{e}_i$ is the quadrilateral value.

The detailed steps of the methodology used to establish the relationship between the hexagonal DGGS and RS images are described in Algorithm 1.

Table

Algorithm 1 Hexagonal pixel modeling
Input: Hexagonal DGGS encoded array:$\mathrm{HexCodeArr}[\mathrm{cod}{e}_1,\mathrm{cod}{e}_2\dots ,\mathrm{cod}{e}_m]$
Output: Hexagonal pixel value $A$
1:	$\mathit{f}\mathrm{or}(i=1\mathrm{to}m)\mathbf{do}$
2:	$\mathrm{HexVertArr}=\mathrm{GetHexVertGeo}(\mathrm{cod}{e}_i)$/* Get the geographic coordinates of the six vertices of a hexagon*/
3:	$\mathit{f}\mathrm{or}(j=1\mathrm{to}6)\mathbf{do}$
4:	$\mathrm{QuadVertArr}=\mathrm{GetQuadVertGeo}(\mathrm{HexVertArr}[j])$/* Calculate the quadrilateral where the vertices of the hexagon are located*/
5:	$\mathit{i}f(\mathrm{ComputeInter}\sec (\mathrm{HexVertArr},\mathrm{quadVetrArr}))==0\mathbf{do}$/*No intersection*/
6:	$A=\mathrm{Valu}{e}_R$
7:	$\mathit{e}\mathrm{lse}\mathbf{do}$
8:	$\mathrm{CalPointInpolygon}()$ /* The points where the hexagon and quadrilateral contain each other, such as $q_1,q_2,h_5$ in Figure 3*/
9:	$\mathrm{CalWeight}()/\ast$Calculate the area of intersecting elements and determine the weight*/
10:	$\mathit{e}\mathrm{nd}\mathrm{if}$
11:	$\mathit{e}\mathrm{nd}\mathrm{for}$
12:	$A$ is calculated using (2)
13:	$\mathit{e}\mathrm{nd}\mathrm{for}$
14:	return $A$

3.1.3. Equal-area and unique mapping between hexagon and rectangle

A high percentage of RS images were arranged in the longitudinal and latitudinal directions, with the centers of the horizontal and vertical pixels parallel to the coordinate axis. The position of the hexagon at any point can be described by skewed coordinate systems (Luczak and Rosenfeld 1976), the symmetric coordinate system (Her 1992; Her 1995; Nagy and Abuhmaidan 2019), and an orthogonal coordinate system. The angles between the skewed coordinate axes are ${\displaystyle \frac{\pi }{3}}$or ${\displaystyle \frac{2}{3}\pi }$, but both have the same essence in 2-D space. We assumed that the two-unit vectors,$v_1$and$v_2$, were linearly independent. The collection of grid center points $L$can be denoted as: (3) $L=\{m{v}_1+n{v}_{2}m,n\in Z\}$(3) Unlike skewed coordinate systems, it is necessary to consider a staggered arrangement of odd and even rows when using an orthogonal coordinate system for indexing. Furthermore, the ratio of the vertical to horizontal distances between adjacent cells is $\sqrt{3}:2$; it is necessary to establish the base vector size in a specific manner to ensure the indexing of any point using integers.

We assumed that the horizontal distance between centers in Figure 4 is denoted as$l$. The base vectors in the horizontal and vertical directions are respectively expressed as$b_1={\left(\begin{array}{cc}{\displaystyle \frac{1}{2}l}& 0\end{array}\right)}^T$ and $b_2={\left(\begin{array}{cc}0& {\displaystyle \frac{\sqrt{3}}{2}l}\end{array}\right)}^T$. Then $L$ is denoted as: (4) $L=m{b}_1+n{b}_2\{\begin{array}{c}m=2k,n=2h(n\text{\%}2=0)\\ m=2k+1,n=2h+1(n\text{\%}2=1)\end{array}k,h\in Z.$(4) Each grid center can be indexed using a 2-D matrix. We can associate the hexagons with quadrilaterals, where the edges of the quadrilaterals are parallel to the directions of the vectors. With this mapping scheme, a unique and equal-area correspondence relationship between them can be established, $S_h=S_q$, where $S_h$ and $S_q$ are the areas of the hexagon and quadrilateral, respectively.

Figure 4. Correspondences between hexagon and quadrilateral. (a) Relationship I. (b) Relationship II. (c) Relationship III.

The three corresponding relationships in Figure 4 indicate that hexagons and rhombuses can be bonded together through relationships I and II. However, we focus on Relationship III because this research is concerned with the organization of RS images for hexagonal pixels. As shown in Figure 5(a), the corner of the rectangle is exactly the center of a hexagon, which can be chosen to represent the rectangular pixels, corresponding to ‘Pixel in corner’ (Baumann 2021). We can alternatively offset the rectangle, so that the hexagonal center is precisely where the center of the rectangular pixel is located (Figure 5(b)); thus, we can represent rectangular pixels by the center point, corresponding to ‘Pixel in center’.

Figure 5. Equal-area correspondence between rectangle and hexagon. (a) Rectangular corner points coincide with hexagonal centers. (b) Rectangular centers coincide with hexagonal centers.

The mapping scheme shown in Figure 5 provides a mathematical foundation for storing spatial projections, coordinate transformations, and position information for hexagonal pixel data using open standard formats.

3.2. Organization of global hexagonal pixel RS images

DGGSs cover the Earth’s surface uniformly with equal-area cells at multiple resolutions. It is generally constructed as a regular grid on the surface of a polyhedron and then projected onto the sphere using an appropriately inverse projection to ensure that cells have an equal-area at the cost of cell shape deformation (Zhou et al. 2020). However, rasters in GIS or RS are primarily organized in a planar rectangular grid with constant distances between adjacent rows and columns.

Developing a direct geometric mapping link between spherical hexagons and rectangles remains a challenge. One approach is to map the hexagon onto the icosahedral surface using a forward polyhedron projection and then flatten it onto a 2-D plane. The equal-area and one-to-one mapping model presented in Section 3.1.3 can then be utilized to enable the conversion from a hexagon to a rectangle (Figure 6).

Figure 6. A unique and equal-area mapping from icosahedral hexagons to rectangles.

3.2.1. Construction of the coordinate systems

Three coordinate systems were established (Figure 7) to illustrate the algorithmic procedure for encoding hexagonal RS data into an open standard format. The specific functions are as follows:

1. The inverse equal-area projection coordinate system,$\mathrm{oxy}$, for which the directions of the coordinate axes are consistent with those of the icosahedral triangular surfaces, allows for the transformation between the plane point and sphere.

2. The 2-D integer coordinate system $\mathrm{oij}$ is defined on a rhombic surface formed by a combination of two triangular surfaces. The cell is identified as $H[d_i,{(i,j)}_n]$, where $d_i$is the index of the rhombus face, $(i,j)$ is the encoding of the hexagonal center, and $n$ is the grid level,$i,j\in Z^{+}$. The encoding is consistent with the idea of using integer indices in regular rectangular grids, which enables efficient neighborhood processing and easy storage.

3. The rhombic Cartesian coordinate system, $\mathrm{OXY}$, which serves two functions. First, it allows for a transitional conversion between the inverse equal-area projection coordinates and 2-D integer coordinates. Second, as the rhombus is formed by a combination of two triangular surfaces, there may be inconsistencies when the study area is located on two different triangular surfaces. The coordinate system also allows the unification of coordinates and facilitates processing.

Figure 7. Illustration of three coordinate systems in an icosahedral rhombus with hexagonal grids.

3.2.2. Hexagonal pixel RS image stored in GeoTIFF

The process of chunking spherical hexagonal RS images into the GeoTIFF is illustrated in Figure 8. First, we projected the geographic coordinates$(\mathrm{la}{t}_i,\mathrm{lo}{n}_i)$of the spherical hexagonal points to the polyhedral surface using the ISEA projection to obtain the inverse equal-area projection coordinate $(x_i,y_i)$of each cell. The rhombus Cartesian coordinate $(X_i,Y_i)$was then calculated based on the $(x_i,y_i)$.

Figure 8. Process of encoding hexagonal data across the rhombic face into GeoTIFF.

Second, we determined the 2-D integer coordinates$(i,j)$ of the geographic boundaries of the specific target object. We assumed the coordinates of $A,C$were $A(i_A,j_A)$ and$C(i_C,j_C)$, respectively. The search range of the matrix for the boundary is$A(0,0),B(0,j),C(i,j),D(i,0)$, where $i=i_C-i_A,j=j_C-j_A$. If the target object extends across multiple rhombus faces, then the searching matrix, $A_k(0,0),B_k(i,0),C_k(i,j),D_k(0,j)k\in [1,10]$, must be calculated separately for each rhombus face. Furthermore, the actual code $(i,j)$ must be adjusted by adding an offset,$A(i_A,j_A)$.

Third, we divided the geographical scope according to the icosahedral rhombus faces. It is expected to traverse $(i,j)$ of each rhombic face separately to determine the rhombic Cartesian coordinates of the smallest bounding $\mathrm{abcd}$:$a(X_1,Y_1),b(X_2,Y_1),c(X_2,Y_2),d(X_1,Y_2).$ Similarly, if the target object extends across multiple rhombus faces, the$a_k(X_1,Y_1),b_k(X_2,Y_1),c_k(X_2,Y_2),d_k(X_1,Y_2),$ for the corresponding rhombus must be calculated.

Finally, the parameters were determined based on the basis vector, $b_1,b_2$, including the image width, image height, and resolution. We can calculate the row and column indices $(I,J)$ of the hexagonal center in the image coordinate system using equations (5). The hexagonal values are then stored in the associated rectangular pixel, along with the projection, coordinate datum and other metadata. (5) $\{\begin{array}{c}I=\lfloor {\displaystyle \frac{\mathrm{Trasform}[5]\times (Y-\mathrm{Transform}[0])-\mathrm{Transform}[2]\times (X-\mathrm{Transform}[3])}{\mathrm{Transform}[1]\times \mathrm{Transform}[5]-\mathrm{Transform}[2]\times \mathrm{Transform}[4]}+0.5}\rfloor \\ J=\lfloor {\displaystyle \frac{\mathrm{Trasform}[1]\times (Y-\mathrm{Transform}[3])-\mathrm{Transform}[4]\times (X-\mathrm{Transform}[0])}{\mathrm{Transform}[1]\times \mathrm{Transform}[5]-\mathrm{Transform}[2]\times \mathrm{Transform}[4]}+0.5}\rfloor \end{array}$(5) Where, $(\mathrm{Transform}[0],\mathrm{Transform}[3])$ is the coordinate of the upper-left point of the image in the $\mathrm{OXY}$ coordinate system,$(\mathrm{Transform}[1],\mathrm{Transform}[5])$ is the image resolution in the$\mathrm{IJ}$direction in $\mathrm{OXY}$, $(\mathrm{Transform}[2],\mathrm{Transform}[4])$ is the angle of rotation of the image, and is the rounding operation.

3.2.3. Multiresolution hexagonal RS data generation

According to the algorithm described above, the connectivity between hexagonal pixels can be organized by regularly arranged rectangles. GeoTIFF can be used in open standard format for storing, transmitting, and visualizing hexagonal data. However, the smallest basic unit remains as a single resolution rectangular structure. Users may still need to implement multiresolution analysis or benefit from hexagons in practical applications. Visualizing or processing all data on a global scale at their original resolution can easily lead to extremely slow computing and consume considerable computational resources. Appropriate multiscale representations can be selected to counteract the adverse impact of increased data volume on the transmission time. A commonly used approach is to parse rectangular pixel images and transform them into refined spherical hexagons at the corresponding level to generate low-resolution grids from high-resolution grids through aggregation, as illustrated in Figure 9. Although the described approach enables the conversion of rectangular pixels into hexagons, it requires a large number of floating-point calculations during the first step.

Figure 9. Multiresolution aggregation of hexagonal data.

A hexagonal aggregation scheme was developed based on the unique mapping model as the aperture-4 subdivision is consistent with the conventional pyramidal structure. The first step is to determine the parent–child affiliation of cells. We refer to the set of children hexagons selected for generation from the seven children hexagons of the parent cell as the child (Figure 10(a)); each refined hexagon at level $n$ has exactly one parent at level $n-1$. As mentioned in Section 3.2.2, hexagons cannot only be arranged using regular rectangles, but various operations can also be performed using optimized matrix calculations on the arrays in memory, including convolution (Onyema et al. 2023), down-sampling, and filtering. Additionally, the unique mapping established allows for the equivalent processing of hexagons while operating on rectangular pixels. In summary, the aggregation methodology can be described as follows: First, we extracted the odd rows and even columns (Figure 10 (b),$A_0$), where all of the reserved ones are the central child of their parent. Then, retrieve six adjacent cells$(A_1,A_2,\dots A_6)$ and determine the value of the parent by weighting the coverage area of the parent and the child. Finally, the cells were restored into sphere. Here, the process of producing multiresolution data, which we refer to as the Hexagonal Multiresolution Pyramid (Hartman and Tanimoto 1984), is similar to the conventional quadrilateral image pyramid. With the help of this scheme, it can equivalently simulate hexagonal LOD by manipulating rectangles, avoiding complex multidimensional aggregation operations. Spherical hexagonal data can be managed using rhombic faces as units shown in Figure 10 (c).

Figure 10. (a) Parent-child affiliation of cells. (b) Hexagonal down-sampling structure. (c) The rhombic surface with aperture-4 subdivided into three consecutive resolutions.

The specific algorithm for computing the spherical grids at level $n=N_{\max }-1$ from level $n=N_{\max }$ is as follows:

Step 1: Take odd rows and even columns of the image as integer coordinate indices$(I,J)$, using equations (6) to calculate the rhombic Cartesian coordinates of the center $(X,Y)$: (6) $\{\begin{array}{c}X=\mathrm{Transform}[0]+J\times \mathrm{Transform}[1]+I\times \mathrm{Transform}[2]\\ Y=\mathrm{Transform}[3]+J\times \mathrm{Transform}[4]+I\times \mathrm{Transform}[5]\end{array}$(6) Step 2: Convert the rhombic Cartesian coordinate $(X,Y)$ to a 2-D integer coordinate,$(i,j{)}_{N_{\max }}$, for the corresponding cell. Determine the coordinates of the parent based on the parent–child relationship and equations (7): (7) $\{\begin{array}{c}{i}_n={\displaystyle \frac{i_{N_{\max }}}{2^{(N_{\max }-n)}}}\\ j_n={\displaystyle \frac{j_{N_{\max }}}{2^{(N_{\max }-n)}}}\end{array}$(7) Step 3: Calculate the projected coordinates, $v_i(i=0,1,2\dots 6)$, of the six vertices based on the projected center coordinates,$(x,y{)}_n$, of the parent. This method depends on the cell position. If the cell lies entirely within a triangular face, it can be calculated directly. However, if the cell falls on the boundary (Figure 9$H_1$), it must be assigned to the corresponding rhombus face before calculation.

Step 4: Determine the value of the parent according to (8). An inverse projection was used to convert the planar grid into a spherical surface, and attribute values of the cells were recorded. (8) $A={\displaystyle \frac{1}{4}{A}_0+{\displaystyle \frac{1}{8}{\sum }_{i=1}^6{A}_i}}$(8) Where $A_0$ is the central child and $A_i$ is neighbor child.

More importantly, the general procedure of such a scheme is theoretically scalable to any DGGS implementation as long as two functions are provided: methods for converting geographic coordinates into grid cells and extracting location information for each vertex.

4. Experiments and analysis

Three experiments were designed to verify the correctness of the mathematical model (Section 3.1.2) and advancement of the storage scheme. Experiment 1 involved hexagonal pixel modeling and precision evaluation. Experiment 2 involved the organization of the hexagonal pixel RS data. To verify the correctness of the mapping model in Section 3.1.3, as well as the feasibility of the multiresolution pyramid, the experiment was divided into two parts: (1) encoding the hexagonal pixel RS data into an open standard format and (2) generating multiresolution hexagonal RS grids from GeoTIFF. Experiment 3 compared the SMOS data organization schemes and visualization. It is used as a comparative object to demonstrate that the solution has more flexibility in sharing and processing SMOS data.

4.1. Experiment 1: hexagonal pixel modeling and accuracy evaluation

4.1.1. Hexagonal pixel modeling

Google’s global data were sampled using an equidistant quadrilateral latitude and longitude grid with a resolution of approximately 1 km, a coordinate reference of WGS84, an image dimension of 21,600$\times$43,200, and a total of 3 bands. However, this geodata framework cannot provide equal-area pixels, which leads to an uneven shape for the global pixels (Kelly and Šavrič 2021; Figure 11(a)). An ISEA4H system was used to organize the global data; the 13th level (average area of 0.76006 km²) was selected according to Table 1, which was closest to the original data. Figure 11(b) shows the results for the Antarctic region.

Figure 11. Diagram of the Antarctic region. (a) The Antarctic region under a latitude and longitude grid. (b) The Antarctic region under a hexagonal grid.

The mean and standard deviation of the compactness (Kimerling et al. 1999) of the cells at different latitudes were determined to assess the uniformity of the hexagon on a global scale. The results are presented in Table 2.

Table 2. Mean and standard deviation of hexagonal compactness at different latitudes.

Grid level	Latitude range
	$0-{30}^{\circ }\mathrm{N}$		${30}^{\circ }\mathrm{N}-{60}^{\circ }\mathrm{N}$		${60}^{\circ }\mathrm{N}-{90}^{\circ }\mathrm{N}$
	Mean	Standard deviation	Mean	Standard deviation	Mean	Standard deviation
6	0.89612	0.00685	0.89723	0.00654	0.89523	0.00664
7	0.89610	0.00692	0.89721	0.00661	0.89520	0.00669
8	0.89610	0.00694	0.89719	0.00664	0.89521	0.00671
9	0.89610	0.00694	0.89719	0.00665	0.89521	0.00671
10	0.89609	0.00695	0.89719	0.00665	0.89521	0.00671

As shown in Table 2, the mean compactness in regions with different latitudes was slightly different at the same level. DGGS uses a hexagon as the basic unit to realize more uniform spatial sampling. In contrast, the latitude-longitude grid can be distorted and deformed at high latitudes, leading to an uneven data distribution. Under the same settings, the standard deviation of cells in the mid-latitude zones was the lowest, demonstrating that these regions had more consistent hexagonal cell variation than the others. The results of the compactness calculations revealed that the geometric characteristics of the cells gradually stabilized as the level increased.

4.1.2. Precision evaluation

In this section, we present a quantitative analysis to demonstrate the advantages of hexagonal pixel modeling. From the perspective of data analysis, information entropy was selected as the basic evaluation strategy for quantifying the information content. The average gradient was selected to measure the clarity of the texture structure. The root mean square error (RMSE) was used to reflect the closeness of the original data.

Entropy reflects the richness of the image information. Comparing the entropy of the modeling results with that of the interpolated sampled images can reflect the degree of information retained by different sampling strategies (Ma et al. 2021). The higher the entropy, the higher the average information content and the higher the image quality, which is calculated as follows: (9) $H(x)=-{\sum }_{i=0}^{N}P(x_i)\mathrm{lo}{\mathrm{g}}_{2}P(x_i)$(9) where $P(x_i)$is the probability of occurrence of pixel grayscale and$N$is the pixel value range.

The average gradient can be used to measure the sharpness of an image and reflect the textural variation of tiny features. The larger the average gradient, the clearer the image. The average gradient formula for the hexagon is presented according to the quadrilateral image calculation strategy: (10) $\mathrm{AverageGrad}={\displaystyle \frac{1}{M\times N}{\sum }_{i=1}^M{\sum }_{j=1}^N\sqrt[3]{{\displaystyle \frac{{\left({\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }x}}\right)}^2+{\left({\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }y}}\right)}^2+{\left({\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }z}}\right)}^2}{3}}}}$(10) where $M\times N$ denotes the total number of pixels, and ${\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }x},{\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }y},{\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }z}}}}$ are derivatives of the three axes in the hexagonal direction.

The RMSE was recommended to indicate the quality of hexagonal images (Ma et al. 2021), calculated as follows: (10) $\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}{\sum }_{i=1}^N{(H_i-R_i)}^2}}$(10) where $H_i$ denotes the hexagonal value, $R_i$ is the rectangular value, and $N$ is the total number of pixels.

Two types of RS images of different resolutions were used. (1) A single-band LandSat8 with a resolution of 30 m released by NASA. (2) Sentinel-2 RGB multiband images were released by the ESA with a resolution of 10 m. Both are UTM projections with a total of 1,000,000 sampling points. The results are presented in Table 3.

Table 3. Results of evaluation indicators from different kinds of images.

Experimental data	Modeling methods	Entropy	Average gradient	RMSE
Sentinel-2	Nearest-neighbor interpolation	7.3709	5.2867	6.5180
	Bilinear interpolation	7.3705	5.7341	6.5846
	Proposed method	7.3894	16.8938	6.2125
Landsat8	Nearest-neighbor interpolation	5.7617	4.2581	5.3302
	Bilinear interpolation	5.7383	4.2326	5.2255
	Proposed method	5.8237	12.7008	4.0425

Based on Table 3, the average gradient of the RS image obtained by the proposed method was approximately three times that obtained by the interpolation methods; there were more textural details and structural features when using the proposed method. Figures 12 and 13 show the results of sampling using different methods at the same level; here, more detailed information, particularly in areas such as rivers, was preserved. The smallest RMSE calculation result obtained in this study indicates that the modeling had the smallest error and was closer to the original data. The proposed method preserved more image information and had a higher clarity. The reason for this can be interpreted as follows: the hexagonal pixel value was calculated based on the distance of a set of neighboring pixels using the interpolation method, where the closer the pixels, the greater their influence on the interpolation result. We chose a more reasonable method in which the attribute value was calculated based on the percentage of area overlap between the quadrilateral and hexagonal pixels, thereby solving the problem of value inaccuracy in the hexagonal pixels. Furthermore, calculating on the sphere avoids the distortion error caused by the projection and preserves their shapes correctly, which leads to a higher modeling accuracy.

Figure 12. Results of the hexagonal pixel modeling of the Sentinel-2. (a) Proposed method. (b) Bilinear interpolation. (c) Nearest-neighbor interpolation.

Figure 13. Results of the hexagonal pixel modeling of the Landsat8. (a) Proposed method. (b) Bilinear interpolation. (c) Nearest-neighbor interpolation.

4.2. Experiment 2: organization of RS images with hexagonal pixels

4.2.1. Storage of hexagonal pixel RS images

This section describes the storage of hexagonal pixel RS data using a standard file format, which simplifies the integration with other map software. We divided global data into 10 rhombuses of equal size. The pixel values and related information were then encapsulated into GeoTIFF. The results of the spherical hexagonal grid storage are shown in Figure 14.

Figure 14. GeoTIFF image corresponding to global hexagonal pixel RS data.

As shown in Figure 14, the global hexagonal pixel RS data were divided into 10 blocks, each of which corresponded to a diamond section on the sphere. Users can utilize GeoTIFF to preview data before obtaining a specific area. In addition, they can combine hexagonal quadtrees and diamond-block images while employing web technologies for direct tile-level data transmission and sharing. This allows rapid data transfer and scheduling for a specific area. If we need to acquire hexagonal data for the entire territory of China, we only need to parse block2 to complete the process quickly, rather than loading global grid data. Some regions are divided into multiple blocks, depending on the orientation scheme between the polyhedron and sphere. The orientation parameters of Zhou et al. (2020) can be referenced if the area of interest must be feasibly stored in one diamond.

In fact, the storage scheme proposed in this paper is compatible with any open standard format, not limited to GeoTIFF, provided that the correspondence between hexagons and rectangles in Section 3.1.3 is established. We can freely choose different files according to our requirements and extend following the storage standards. Figure 15 shows the results of storing global hexagonal grid data using HDF5. The hexagonal pixel values are stored as matrices in different file formats, with the main difference being the metadata management. GeoTIFF adopts a series of Keys to access the information, but HDF5 can add metadata directly to describe information in the file, group or dataset.

Figure 15. Hexagonal pixel images were stored using HDF5.

Owing to the arrangement of the hexagon, images on rectangular grids are also presented to users in an alternating odd-even row structure, which is different from the traditional arrangement of square-pixel images. We refer to this as the simulated hexagonal dataset. The open file format can be easily interfaced with common GIS applications, libraries, and tools to manage attribute values, spatial projections, and coordinate information of hexagonal images in a standardized manner.

4.2.2. Multiresolution hexagonal RS data

In this section, an aggregation scheme is used to generate hexagonal grids of different resolutions. We began with the data (Figure 14) at resolution of n, after which the resulting image array was decimated to leave only the values in the cells at a resolution of n-1. This process was repeated until the desired data resolution was achieved. Subsequently, spherical hexagonal data of different resolutions can be generated using these parameters or metadata as shown in Figure 16.

Figure 16. RS data for spherical hexagonal pixels at different resolutions.

During our experiment, there were several beneficial characteristics of the scheme.

1. The aperture-4 hexagon subdivision and the unique mapping scheme can be combined to process hexagonal pixels in the form of arrays in memory, decreasing the complex floating-point calculation, which improves the computational efficiency.

2. The parent–child addressing relationship constructed based on the hexagonal DGGS can be calculated quickly using binary operations, as illustrated in equation (6(6) $\{\begin{array}{c}X=\mathrm{Transform}[0]+J\times \mathrm{Transform}[1]+I\times \mathrm{Transform}[2]\\ Y=\mathrm{Transform}[3]+J\times \mathrm{Transform}[4]+I\times \mathrm{Transform}[5]\end{array}$(6) ), which reduces the complexity of the parent–child query.

3. The encapsulation of recovery parameters reduces the computational resources required for the conversion of coordinate systems during data processing.

4.3. Experiment 3: comparison of SMOS schemes

4.3.1. Organization of SMOS

The SMOS mission, scheduled for launch in early 2009, was part of the ESA Living Planet Programmer (Martin Suess 2004). The SMOS data were sampled using the ISEA4H DGGS, which presents uniform adjacency with a mean intercell distance of approximately 15.0 km. The scientific data for SMOS L1 and L2 are available at (SMOS Online Dissemination 2022). The true data were restored according to the geographic coordinates and salinity values extracted from four consecutive days of ocean salinity (OS) data stitching in July 2022. We added global continental coastline lines to verify the accuracy of the modeling results for the geographical locations shown in Figure 17.

Figure 17. Global hexagon pixel OS data with continental coastline lines.

The SNAP software provided by the ESA allows for a quick preview of the hexagonal data, which simulates an approximate hexagon through multiple square pixels. Although this facilitates the use of existing algorithms based on the rectangular pixels, it significantly increases the volume of data, and the precision of the representation is inaccurate Furthermore, it is important to note that the number of square pixels forming a single hexagon is not constant (averaging approximately 40 pixels), as depicted in Figure 18(a). Consequently, this results in an average size of 5.63 km² for each square pixel. To make a more reasonable comparison, the data were resampled to the 12th level (with an average area of approximately 3.04 km²) and stored in the GeoTIFF format, as illustrated in Figure 18(b).

Figure 18. Image details of the two schemes for the OS data. (a) SNAP. (b) Proposed method.

Both schemes can simulate hexagonal OS data. Nevertheless, it is worth noting that the templates employed for simulating the hexagon were not immutable. Consequently, accurately determining the position of a simulated hexagon within a rectangular pixel image has become a challenging endeavor, thus hindering subsequent computational tasks. In contrast, the morphology of the pixels in this study remained constant, thereby facilitating the formulation of fixed convolutional operation templates for image processing.

We ensure a precise correspondence in position and geometric equivalence in the area when transforming hexagons into rectangles, which has significant practical implications for data processing compared to ESA. Precise correspondence in the realm of physical space preserves the inherent geometric attributes of the hexagonal staggered arrangement. This enables the prompt identification of an adjacent cell within the spatial domain. Equivalency in terms of area allows for direct statistical analysis based on the number of rectangular pixels and attribute values, without considering the influence of cell sizes. This can be logically considered to process hexagonal pixel image. Furthermore, we also have an advantage in processing efficiency, as our approach handles a single rectangular pixel equivalently as a hexagon, whereas current SMOS scheme requires processing multiple square pixels. Data opening means that the data is published publicly and available for open access, which serves as the premise and foundation of data sharing. We adeptly organize and manage hexagon OS data in an open file format, thus facilitating its further sharing.

The ESA offers data collection using a pure aperture. Resampling algorithms are likely be adopted to produce multiresolution hexagonal data. However, this destroys the correspondence between the hexagon and multiple squares, making it difficult to convert square pixels into hexagonal grids. Instead, we can iteratively obtain GeoTIFF images corresponding to different grid levels and then convert them into real hexagons. Figures 19 and 20 show the OS data for the rectangular and hexagonal pixels at three consecutive resolutions in layers 8, 9, and 10. The hexagonal pixels maintained an accurate positional correspondence with the rectangular pixels in the same region without any error. This scheme enhances the flexibility and adaptability of the approach in practical applications, making it a viable option for the organization of SMOS data.

Figure 19. Rectangular pixel OS data at different resolutions. (a) Level 8. (b) Level 9. (c) Level 10.

Figure 20. Hexagonal pixel OS data at different resolutions. (a) Level 8. (b) Level 9. (c) Level 10.

4.3.2. Visualization with ArcGlobe and WorldWind

ArcGlobe is a component of the ArcGIS desktop system's 3-D extension. It offers continuous, multiresolution, and interactive viewing of global geographic data used across many sectors and industries and supports the reading and parsing of multiple file formats. Consequently, we decided to import the simulated hexagonal OS data into the ArcGlobe to provide users with a visualization. WorldWind (World Wind 2023) is a digital Earth model provided by NASA that is completely open-source, accessible, and highly scalable. It also enables engineers to create application tools for viewing, manipulating, and analyzing spatial data. Therefore, we chose WorldWind to load and display the hexagonal data.

One of the strips of data was selected for experimentation and encoded into a GeoTIFF file. However, the GeoTIFF geographic codes are based on reference information from the European Petroleum Survey Group (EPSG) database, which does not currently support the ISEA coordinate system. We can also customize the polyhedral projection metadata in the form of WKT, but the GIS/RS software cannot be recognized directly. We further converted the ISEA projection coordinate system to an EPSG-supported georeferencing system (e.g.4326:WGS84) and provided users with data visualization using ArcGlobe (Figure 21). The scheme presented in this paper is more flexible. Open file formats such as GeoTIFF allow users to parse the stored hexagonal data according to standards. In general application scenarios where there is no need to restore the hexagonal cells to spherical surface, the hexagonal OS dataset can be used as a regular image. In specialized application scenarios, dedicated algorithms can accurately restore it to a spherical surface based on the stored metadata information. Figure 22 shows a visualization of the hexagonal data for the same area using WorldWind.

Figure 21. Visualizing rectangular pixel OS data using ArcGlobe. (a) Global view. (b) Local view.

Figure 22. Visualizing rectangular pixel OS data using WorldWind. (a) Global view. (b) Local view.

5. Conclusions and future work

An aperture-4 hexagonal DGGS was used as the bearer data framework of the RS, through which we established a mathematical model to obtain a high-precision representation and efficient organization of hexagonal pixel RS images. Our main findings were as follows. First, as demonstrated in Experiment 1, the modeling scheme had higher accuracy and retained more structural features. The data in different ranges could be uniformly represented on a sphere. Second, an equal-area and one-to-one mapping relationship established between the hexagon and rectangle allowed for the management of hexagonal pixel images with open standard formats, which can make the data compatible with the current common GIS software. In addition, complex operations on hexagons can be converted into efficient matrix operations. The hexagonal LOD could be equivalently simulated by manipulating the rectangles. Finally, we compared the scheme with SMOS and found that not only can we achieve the simulation of hexagonal pixel images, but the equal-area correspondence established also made it more flexible for acquiring different-resolution data.

This study provides guidance for modeling, organizing, and sharing hexagonal pixel images, which can facilitate broader extensions of DGGS. Future work should focus on algorithmic processing of hexagonal pixel images and discover more information for mining, such as in field of human emotions recognition (Chhabra et al. 2022).

Acknowledgments

We would like to express our great appreciation to the editors and five anonymous reviewers for constructive comments that helped improve the manuscript.

Disclosure statement

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

Data availability

The data in the experiment can be downloaded from: https://smos-diss.eo.esa.int/oads/access/ and https://www.gscloud.cn/.

Additional information

Funding

This work was supported by the Special Science Fund for Innovation Ecosystem Construction of National Supercomputing Center in Zhengzhou [grant number: 201400210100].