1. Introduction
Organizing and referencing geospatial data pose increasingly complex challenges due to the sheer volume and rapid growth of the data collected. For example, the daily influx of highresolution satellite imagery alone amounts to terabytes of data. It is important to store this data in a format that allows easy access, referencing, sharing, and analysis without frequent reprojections to maintain accuracy.
Addressing these challenges requires the development of a spatial reference frame capable of fusing data from diverse sources into a global mosaic at multiple resolutions. Discrete Global Grid Systems (DGGSs) have emerged as a promising class of such reference frames. They use hierarchical tessellation to partition and address the entire planet without gaps or overlaps. The development of DGGS began in the mid20th century [
1], but only became popular at the turn of the century [
2,
3,
4]. A prominent subclass of DGGSs are Geodesic DGGSs (GDGGSs) [
4], which project the surface of the planet onto the faces of regular or semiregular circumscribed polyhedra. Commonly used base polyhedra are the five Platonic solids, in particular the hexahedron (cube), while the most common semiregular polyhedron is the truncated icosahedron. The appearance of unfolded regular polyhedra, the number and shape of faces, and the way they partition the Earth's surface can be seen in
Figure 1.
While increasing the number of faces in a polyhedron improves its approximation of a spherical surface, it introduces complications when merging adjacent partitions due to separate coordinate systems and data sets. Moreover, memory and CPU consumption increase with the number of partitions in systems handling streaming geospatial data, like outofcore planetscale terrain rendering applications [
5]. Therefore, minimizing the number of partitions becomes desirable. This paper explores the use of orthogonal cylindrical surfaces, also known as the YinYang grid [
6], as projection planes to reduce the partitions to two. While the YinYang grid has found applications in various fields, its potential as a cartographic projection has not been sufficiently explored. Additionally, this paper investigates the rotation of projection cylinders to reduce distortion of Earth's landmass and minimize disruption of continental plates caused by partition boundaries. The utilization of auxiliary latitudes in reducing distortion during the mapping of the ellipsoid onto the sphere is also considered.
This paper is divided into five sections. After this brief introduction,
Section 2 provides an overview of the historical development and importance of organizing global geospatial data, as well as the emergence and standardization of DGGSs.
Section 3 addresses the projection of geospatial data onto a plane using the proposed method, which encompasses ellipsoid to sphere mapping, sphere to orthogonal cylinders mapping, and the orientation of orthogonal cylinders. The goal is to minimize distortions and the number of partitions. Experimental results and discussion are presented in
Section 4, followed by the conclusion in
Section 5.
2. Related Work
The need to develop a system for organizing global geospatial data is not new. One of the first studies of the feasibility of implementing the Erath Data Base System [
7] was conducted in the early 1970s for the needs of the U.S. Navy. The system was based on the Quadrilateralized Spherical Cube (QSC), one of the first hexahedral projections implemented on digital computers. The proposed system was soon modified [
8], and in the following years it was also used as part of the Cosmic Background Explorer (COBE) project at NASA.
Due to the regular and uniform structure of the grid consisting of square cells, consistency with the Cartesian coordinate system, ease of interpolation and extrapolation, and straightforward visualization, hexahedral projections have gained wide popularity and are used in many different fields.
Hierarchical Equal Area isoLatitude Pixelization (HEALPix) is a class of spherical projections with the property of distributing 12N
^{2} points as uniformly as possible over the surface of the unit sphere [
9]. These hybrid projections combine the Lambert cylindrical equalarea projection for the equatorial region with the interrupted Collignon projection for the polar regions. Of this infinite class of projections, only the projection with 3 base resolution pixel layers between the north and south poles and 4 equatorial base resolution pixels can be rearranged to a hexahedral projection.
Rotated HEALPix (rHEALPix) [
10] is an extension of the HEALPix scheme that introduces rotation capabilities, is better adapted to standards, and inherently combines polar triangles into quadratic partitions. rHEALPix has found wide application in organizing global geospatial data [
10,
11,
12]. All of the previously mentioned projections are equalarea, but with significant angular distortions and even discontinuities. Due to their simpler implementation and relatively good balance between angular and areal distortion, many hexahedral projections also find application in computer graphics [
13] and, in particular, in the planetscale terrain visualization [
14]. Among the best known are: Adjusted Spherical Cube (ASC) [
15], Continuous Cube Mapping (CCM) [
16], and Cartesian Spherical Cube (CSC) [
17].
The expansion of the system for organizing and referencing global geospatial data occurs at the turn of the century, when the first classifications appear and the term Discrete Global Grid System (DGGS) is introduced for a very significant class of such systems. The importance of DGGSs is also reflected in the fact that the Open Geospatial Consortium (OGC) established the DGGS Standard and Domain Working Groups to support the standardization of these systems. In 2017, the OGC published the first version of the DGGS Abstract Specification [
18]. The standardization process continued, resulting in a formal specification defined by the ISO 191701:2021 standard [
19] and a revised version of the OGC Abstract Specification [
20]. The popularity of DGGS has also increased due to numerous open source implementations. [
21]
A recent trend in geospatial data processing involves the implementation of datacubes based on DGGSs [
22], enabling efficient management of big data workflows. DGGSs, as a standardized representation of the Earth, provide the foundational platform for Digital Earth [
23]. Digital Earth is a concept that aims to create an interactive digital replica of the entire planet, fostering a shared understanding of the relationships between the physical and natural environment and society [
24].
The design of GDGGSs is characterized by five fundamental elements [
4]:
A regular base polyhedron;
The orientation of the base polyhedron with respect to the planet;
A hierarchical spatial partitioning of the polyhedron faces;
The mapping of a spherical or ellipsoidal surface to polyhedral faces and vice versa;
Methods for indexing and addressing cells.
In the next section, we propose improvements to three of these properties of GDGGSs, aiming to reduce the number of partitions while minimizing distortion effects. Specifically, we replace the regular polyhedron with orthogonal cylinders whose orientation is defined to minimize landmass distortion. We use an equidistant cylindrical projection on two orthogonal cylinders, called the YinYang grid [
6], to map the Earth's surface. Additionally, we consider the use of approximated auxiliary latitudes in mapping the ellipsoid onto the sphere.
3. Method Description
Organizing geospatial data of planet Earth is a significant challenge, particularly with the recent influx of large volumes of data from diverse sensors that require integration into a cohesive mosaic while ensuring accessibility and interoperability. This data organization should be efficient in terms of:
Storage–utilizing a compact distributed approach.
Acquisition–enabling simple addressing and supporting spatial and temporal locality.
Analysis–ensuring data is in a suitable form for processing, preferably without the need for reprojection during usage and with minimal loss of precision in transformations.
Visualization–storing data in a format suitable for display.
Meeting all these requirements simultaneously is not easy, given the diverse applications for geospatial data. Many applications focus solely on the Earth's surface, with data recorded as twodimensional raster layers, thereby determining the shape of the referencing system. However, projecting the spheroidal surface of the planet onto a plane has been a longstanding challenge for cartographers. No transformation fully preserves all the properties of surface entities, ensuring that they retain their shape, proportionality, and continuity. Conformal projections preserve shape but distort area significantly, while equalarea projections preserve area but distort shape considerably. A projection that is both conformal and equalarea does not exist. Preserving one property more effectively comes at the expense of the other. Furthermore, representing the entire planet's surface on a single plane without singularities and discontinuities is not possible.
The goal of the reference frame is to establish a uniform tessellation of the planet's surface with no gaps and a unique cell addressing system. Due to the planar organization of data, it becomes necessary to partition the surface into multiple sections. Increasing the number of partitions reduces distortion but complicates data manipulation when merging two or more partitions, as each partition employs its own local coordinate system. These partitions are further divided into sections, which consist of blocks of data suitable for retrieval and processing. Subsequently, the sections are subdivided into smaller units known as cells, which represent the smallest addressable units in the system.
Figure 2 illustrates the process of transforming the planet's surface into an addressable system of cells, using the example of a subdivision into two partitions based on orthogonal equidistant cylindrical projections.
3.1. Mapping an Ellipsoid to a Sphere
The ellipsoid is the most commonly used approximation for the shape of the planet Earth today. Due to the widespread utilization of the Global Positioning System (GPS) and the abundance of data collected using this reference frame, the WGS 84 ellipsoid [
25] serves as the primary model for creating global datasets. To ensure interoperability, the NGA (National GeospatialIntelligence Agency) closely aligns the WGS 84 reference frame with other standards, particularly the International Terrestrial Reference Frame (ITRF) [
26]. Consequently, the latest official revision of the WGS 84 reference frame (G2139) remains consistent with the IGb14 realization of the ITRF2014 [
27].
In geodetic computations, the ellipsoid model is often substituted with the spherical model due to its higher symmetry and simpler calculation of many geodetic formulas. Since the ellipsoidal model deviates slightly from a perfect sphere in the case of the Earth, the spherical formulas can be applied to the ellipsoid by replacing the geodetic latitude by one of the "auxiliary latitudes". Introducing a mapping from an ellipsoid to a sphere introduces an additional distortion that varies depending on the applied auxiliary latitude. There are six auxiliary latitudes: geocentric, conformal, authalic, parametric, rectifying, and isometric. O. Adams systematically described these auxiliary latitudes and derived all the corresponding formulas in 1921 [
28], but they gained popularity later after the publication of Snyder's working manual [
19].
The basic latitude used in global datasets and position determination based on global navigation is
geodetic latitude (θ). It represents the angle between the equatorial plane and the surface normal at a point on the ellipsoid. Calculating geocentric latitude (φ), which represents the angle between the equatorial plane and the radius vector, is relatively simple compared to other auxiliary latitudes. Equation (1) can be used to calculate the geocentric latitude based on the geodetic latitude, where
e denotes the eccentricity of the ellipsoid.
Two auxiliary latitudes are of significant importance in addressing specific types of distortion. The application of the conformal latitude (χ) results in a conformal mapping of an ellipsoid onto a sphere, effectively eliminating angular distortion. On the other hand, the use of the authalic latitude (β) achieves an equalarea mapping, eliminating areal distortion. The conformal latitude can be computed from the geodetic latitude using equation (2).
Computing the geodetic from the conformal latitude, i.e. the inverse transformation, requires an iterative procedure or series [
28,
29]. Equation (3) is one of the methods for the inverse transformation. By using the first four terms of the sum, the computational error can be kept below 1.0E12. The corresponding values for the coefficients
c_{i} can be found in [
29].
The authalic latitude is calculated from the geodetic latitude using equations (46).
Similar to the conformal latitude, the inverse transformation for the authalic latitude requires an iterative procedure or series.
Figure 3(a) depicts the deviation of the aforementioned auxiliary latitudes from the geodetic latitude. It can be observed that the geocentric and conformal latitudes have very similar deviations from the geodetic latitude.
Figure 3(b) displays the difference between the conformal and geocentric latitude values as a function of geodetic latitude. The maximum deviation occurs at 60° north and south (geodetic) latitude and is approximately 1.4E2° (50.4"). Considering the small deviation from conformal latitude, the relative simplicity of the calculation, and the availability of a straightforward closedform inverse transformation, geocentric latitude can be used to mitigate the angular distortion of mapping an ellipsoid onto a sphere.
To simplify the calculation of authalic latitude and gain closedform inverse transformation, we propose the following approximated formula:
The smallest maximum deviation of the approximation (β') from the authalic latitude (β) for the WGS 84 ellipsoid is obtained with k = 0.666741 and is smaller than 2.16E5° (0.078"), as shown in
Figure 3(c). This represents a 33% improvement over the approximation presented in [
30]. The proposed formula is similar to the geocentric latitude formula and both the forward and inverse transformations have closed forms and can be easily computed. The ratio of the tangents of the geodetic latitude to the approximated authalic latitude for the WGS 84 ellipsoid is approximately 1.004488.
3.2. Mapping a Sphere Onto Orthogonal Cylinders
While a sphere is more suitable for geodesic calculations compared to an ellipsoid, it cannot be flattened into a plane without interruptions. To overcome this limitation, the next step involves projecting the sphere onto another figure that has flat surfaces or can be unrolled into a plane seamlessly. GDGGSs achieve this by utilizing the faces of circumscribed regular or semiregular polyhedra as the projection surfaces for the sphere.
Each face of the polyhedron represents a partition of a specific projection. Platonic solids are commonly used as the base polyhedral [
4] because they possess regularity, with faces of the same shape (triangles, squares, or pentagons), equal areas, and an equal number of neighboring faces. Among the regular polyhedra, the hexahedron (cube) is particularly popular due to its relatively small number of partitions (six) and the square shape of both the partitions and cells.
In addition to regular polyhedra, a commonly used semiregular polyhedron is the truncated icosahedron [
31], which belongs to the group of 14 Archimedean solids. The truncated icosahedron lends itself to a hexagonal cell structure, although its sides are not all identical. It consists of 12 pentagonal and 20 hexagonal sides.
Cylinders are often used as an alternative to flat surfaces for projections. Due to its better adhesion to the spherical surface, the cylinder produces less distortion even with a smaller number of projection surfaces. One of the oldest and simplest cylindrical projections still in use today is the equidistant cylindrical projection. It was invented by Marinus of Tire around 100 AD and, despite its considerable distortion, remains the most commonly used projection for organizing global data. This projection establishes a straightforward relationship between map positions and the corresponding geographic locations.
By combining two orthogonal equidistant cylindrical projections, a projection known as YinYang [
6] is obtained. The YinYang projection, together with its associated grid, has found applications in a variety of fields, such as simulations of geodynamo and mantle convection [
6], visualization of 3D mantle convection [
32], global shallow water models [
33], 3D hydrodynamic simulations of corecollapse supernova evolution [
34], feature extraction from omnidirectional panoramic images [
35], and many others. However, its potential as a cartographic projection has not been extensively explored.
In general, the YinYang approach uses two complementary components, and the mapping need not be based on orthogonal equidistant cylindrical projections. However, for simplicity, it is usually implemented in this way. To refer to the projection more precisely in this paper, we use the term Dual Orthogonal Equidistant Cylindrical (DOEC) projection.
The first partition (P0) generally extends along the equator and is symmetrical about the equator and the prime meridian. It uses polar coordinates identical to global geographic coordinates (ϕ, θ). The second partition (P1) extends along the antimeridian (180
^{th} meridian), is symmetric about it, and includes both poles.
Figure 2(d) shows the two partitions in a rectangular shape that facilitates the determination of their boundaries. In this shape, however, the partitions overlap by about 6.4% of their total area.
Figure 4 shows the partitions without overlaps and their distribution over the globe.
The local coordinate system of partition P0 coincides with the global geographic grid, so no coordinate conversion is required. On the other hand, the local polar coordinates of partition P1 are determined based on either the geographic coordinates or the coordinates of partition P0 using equations (8) and (9).
Because of the orthogonality of the partitions, the formulas (8) and (9) can also be used to convert coordinates from partition P0 to P1 by simply exchanging the arguments. The condition indicating that a point with local polar coordinates (ϕ
_{p}, θ
_{p}) belongs to the current rectangular partition, including overlapping areas, is defined by the logical expression (9). The expression (10) additionally indicates that the point belongs to the overlapping area, where θ
_{q} is the latitude in the complementary partition obtained by equation (8).
3.3. Orientation of Projection Cylinders
The distribution of the distortion depends on the projection applied and is never uniformly distributed over the surface of the partition. Typically, the distortion is minimal in the center of the partition and increases toward the edges and corners, indicating a greater distance of the projection plane from the surface of the sphere.
Tissot's indicatrices [
29] are commonly used to visualize distortions. They are represented as ellipses formed by projecting infinitesimal circles from the surface of the globe onto the projection plane. The size, eccentricity and inclination of these ellipses indicate the type and degree of distortion present.
To better observe the distribution of the deformation parameters, instead of using ellipses that combine multiple deformation parameters graphically, we represent each parameter individually using a red color intensity scale. The values of the distortion parameters are still calculated using the indicatrices. For instance, the angular distortion is determined by calculating the maximum angular deformation ω using equation (12), based on the major (a) and minor (b) semiaxes of the indicatrix. The formulas for calculating the indicatrices and the corresponding deformation parameters can be found in [
29].
Figure 5 illustrates the distribution of angular, areal, and aspect distortion for both DOEC partitions for the basic orientation of the projection cylinders, as described in the previous section. The aspect distortion refers to the ratio of the major and minor semiaxes of the indicatrix (a/b). For DOEC projection, aspect distortion can be interpreted as normalized surface distortion, where the minimum value means no distortion.
By changing the orientation of the cylinders, we can affect the distribution of feature distortion on the planetary surface. Changing the position of the projection planes or circumscribed polyhedron [
36] is a common technique in cartography to achieve certain desired effects. Even in the oldest hexahedral projection [
37], the base cube is rotated 27° about the Earth's axis of rotation to align the westernmost point of the African continent with the edge of the cube. An early discussion of the orientation of the base cube can be found in [
8], but without an indepth examination of the underlying considerations and applications.
The main reasons for the change of the basic orientation, which assumes alignment with the equator, the prime meridian and the poles, can be summarized as follows:
Avoiding fragmentation of target areas: Adjusting the orientation helps prevent splitting local or regional target areas across multiple faces of the polyhedron [
38]. This ensures the integrity of these areas in the projection.
Encompassing an entire continent: Changing the orientation allows an entire continent, such as North America, to be included in a single partition [
8,
31]. This is beneficial for regionally focused mapping and analysis, and is suitable for the polyhedral with fewer and larger faces.
Preventing ruptures in the continental plates after the base polyhedron has unfolded: This is achieved by positioning the vertices of the polyhedron at the oceans, as demonstrated in Fuller's Dymaxion Airocean World Map [
39].
Minimizing landmass distortion: Another important criterion for the orientation is minimizing landmass distortion [
30]. This aims to preserve the accurate representation of land features on the map.
The first two reasons mentioned above are locationspecific and may not be applicable to global systems since they prioritize specific regions. On the other hand, the latter two reasons are more universal and serve as criteria for determining the orientation used in DOEC. However, both criteria cannot be met at the same time. So, we first used an iterative process and a rotation around all three Cartesian axes with a check of angular and areal distortions to determine the optimal orientation considering landmass distortion.
To determine the optimal orientation, a vector map of the world [
40] was rasterized in LatLon WGS 84 (EPSG:4326) projection at a resolution of 4096×2048 pixels (see
Figure 6). The map was then reprojected onthefly into two partitions with progressively varying rotation angles. Nearest neighbor sampling was used in the reprojection since it is fast and clearly delineates the continents. In addition, Antarctica was excluded from the map to focus on more cartographically important areas. The distortion is checked only for the pixels that belong to the landmass.
Considering the rotation angles, denoted as ϕ
_{r} for the vertical axis (longitudinal rotation), θ
_{r} for the horizontal axis (latitudinal rotation), and ρ
_{r} for counterclockwise rotation about the axis perpendicular to the previous two, the minimum distortion of the landmass at one degree resolution was obtained as ϕ
_{r} = 131°, θ
_{r} = 49°, and ρ
_{r} = 20°. To further reduce the clipping of the continental plates, the following corrected rotation angles are proposed: ϕ
_{r} = 125°, θ
_{r} = 50°, and ρ
_{r} = 15°.
Figure 7 illustrates the layout of the partitions based on the proposed rotation angles. The advantages of DOEC projection and the effects of the proposed method on distortion reduction are shown in the next section.
4. Experimental Results and Discussion
The DOEC projection offers several key advantages over circumscribed polyhedra, the typical choice in DGGSs. First, it reduces the number of partitions to only two. This property proves beneficial for outofcore terrain rendering algorithms such as Ellipsoidal Clipmaps [
5]. Minimizing the number of partitions displayed simultaneously reduces memory consumption, since each partition requires corresponding structures for visualization. These structures typically include terrain elevations and highresolution aerial imagery at multiple levels of detail. For outofcore algorithms, the constant updating of these structures with geospatial data places a burden on the central processor and requires access to network resources or slow secondary media to retrieve the data, so the number of partitions directly affects system performance.
In addition, partition connections require special treatment in threedimensional visualization that includes additional testing of conditions, clipping, and fitting. Therefore, minimizing the occurrence of partition connections improves overall performance. Moving the partition boundaries above the water surface further facilitates seamless joining. Consequently, the proposed DOEC projection incorporates rotated projection cylinders to not only reduce distortion, but also minimize clipping of the continental plates by the partition boundaries.
Another advantage of the DOEC projection is its ability to achieve a favorable balance between areal and angular distortions.
Table 1 shows the distortion values for common projections, considering a perfect sphere and without rotations of the projection planes. These values were derived from more than 67 million measurement points evenly distributed over the surface of the partitions.
DOEC is compared to other projections that also use square cells and can be treated as hexahedral projections. These include Quadrilateralized Spherical Cube (QSC), rotated Hierarchical Equal Area isoLatitude Pixelization (rHEALPix), Adjusted Spherical Cube (ASC), Continuous Cube Mapping (CCM), and Cartesian Spherical Cube (CSC).
rHEALPix is a hybrid projection that combines the Lambert cylindrical equalarea projection for the equatorial region with the interrupted Collignon projection for the polar regions. To distinguish between these two regions, they are referred to as rHEALPix
^{E} and rHEALPix
^{P}, respectively, in
Table 1. Both the QSC and rHEALPix
^{P} projections are equalarea projections, but have significant angular distortions and even discontinuities along the diagonals of the partitions. This is due to the fact that their partitions consist of four triangular surfaces. The other three projections, on the other hand, offer a more balanced compromise between surface distortion, angular distortion, and ease of implementation.
Figure 8 provides a visual comparison of the projections listed in
Table 1, offering insights into the distortion effects. Although the rHEALPix
^{P} projection is primarily designed for the polar regions, it is applied to a portion of the equatorial region (
Figure 8(f)) to demonstrate the distortion effects on similar shapes as the other projections. To quantify the combined effect of angular and areal distortions,
Table 1 contains the geometric mean of the average values of angular and normalized areal distortions calculated using equation (13).
In equation (13), the normalized averaged areal distortion (σ
_{ave/min}) is used instead of the simple average (σ
_{ave}) because it better captures the variation in area change and obtains a value of 1 when no distortion is present.
Table 1 shows that the DOEC projection has the lowest values for maximum and average angular distortion, maximum and average aspect distortion, and relatively low normalized averaged area distortion. Consequently, it achieves the lowest GM
_{ωσ} value among all the projections presented.
Table 2 provides insight into the effects of projection cylinder rotation on average distortions for some characteristic values. The values listed in the table were obtained by averaging over 6.7 million measurement points belonging to the landmass. The location of the landmass was determined by reprojecting a rasterized world vector map [
40] using the appropriate rotation of the projection cylinders. Applying the rotations ϕ
_{r} = 125°, θ
_{r} = 50°, and ρ
_{r} = 15° (hereafter abbreviated as R(125,50,15)) reduces the average angular distortion of the landmass by 1.9 times and the normalized average areal distortion by about 6%.
The rotation R(131,49,20) yields the lowest average value of angular distortion, but causes the partition boundary to intersect the southern part of the African continent (
Figure 9), resulting in a reduction in average angular distortion of less than 1%. Considering the optimization of the disruptions in the continental plates, the rotation R(125,50,15) proves to be a superior solution. It effectively reduces the ruptures in the continental plates while providing a notable improvement in the average distortions.
Table 3 presents a comparison of the distortion parameters for different authalic latitudes applied to the entire surface of partition P0 without rotating the projection cylinders. The results show that applying conformal latitude to map an ellipsoid onto a sphere yields identical values for angular and aspect distortion as for an ideal sphere (see first and second rows in
Table 3), indicating that this mapping does not result in any additional angular distortion. The application of the authalic latitude also shows no additional areal distortion. The average angular distortion deviation when using the geocentric latitude is about 0.0026%, while the average normalized areal distortion using approximate authalic latitude is less than 0.00018%. These results indicate that the use of approximated auxiliary latitudes does not introduce significant additional distortion. It is worth noting that the use of approximated authalic latitude even reduces the average angular distortion in this particular case. Its simplicity and favorable performance make it the optimal choice for mapping the ellipsoid onto the sphere.
All previously presented results indicate that the proposed mapping of Earth's geospatial data can be effectively used to mitigate distortions and organize the data into only two partitions.
5. Conclusions
The Dual Orthogonal Equidistant Cylindrical (DOEC) projection offers several advantages over other projections in terms of reducing memory consumption and improving overall performance for worldscale geospatial data visualization, minimizing partition interconnections, and achieving a favorable balance between areal and angular distortions. Compared to other projections such as Quadrilateralized Spherical Cube (QSC), rotated Hierarchical Equal Area isoLatitude Pixelization (rHEALPix), Adjusted Spherical Cube (ASC), Continuous Cube Mapping (CCM), and Cartesian Spherical Cube (CSC), the DOEC projection has lower values for maximum and average angular distortion, maximum and average aspect distortion, and relatively low normalized average area distortion. It also incorporates rotated projection cylinders to minimize landmass distortion and continental plate disruption. The optimal rotation of the projection cylinders reduces the average angular landmass distortion to about 3.6° and the normalized average area distortion to about 1.07, while optimizing continental plate disruptions. Furthermore, by applying the approximated authalic latitude, the DOEC projection preserves the areal distortion while additionally reducing the average angular distortion and is considered the best candidate for mapping an ellipsoid onto a sphere. Overall, the DOEC projection provides improved performance and distortion characteristics, making it a valuable choice for mapping global geospatial data.
Author Contributions
Conceptualization, A.D.; methodology, A.D; validation, A.M.; writing—original draft preparation, A.D; writing—review & editing, A.D. and A.M.; visualization, A.D.; supervision, D.R.; funding acquisition, D.R.
Funding
This research was supported by Ministry of Education, Science and Technological Development of the Republic of Serbia, grant number 4510347/202301/200102.
Conflicts of Interest
The authors declare no conflict of interest.
References
 Vestine, E.H.; Sibley, W.L.; Kern, J.W.; Carlstedt, J.L. Integral and SphericalHarmonic Analyses of the Geomagnetic Field for 1955.0, Part 2. J. Geomagn. Geoelectr. 1963, 15, 73–89. [Google Scholar] [CrossRef]
 Kimerling, A.J.; Sahr, K.; White, D.; Song, L. Comparing Geometrical Properties of Global Grids. Cartogr. Geogr. Inf. Sci. 1999, 26, 271–288. [Google Scholar] [CrossRef]
 Alborzi, H.; Samet, H. Augmenting SAND with a spherical data model. In Proc. First Int. Conf. on Discrete Global Grids; Santa Barbara, CA, USA, 2000; pp. 26–28. [Google Scholar]
 Sahr, K.; White, D.; Kimerling, A.J. Geodesic Discrete Global Grid Systems. Cartogr Geogr Inf Sci 2003, 30, 121–134. [Google Scholar] [CrossRef]
 Dimitrijević, A.; Rančić, D. Ellipsoidal Clipmaps–A planetsized terrain rendering algorithm. Computers & Graphics 2015, 52, 43–61. [Google Scholar] [CrossRef]
 Kageyama, A.; Sato, T. YinYang grid: An overset grid in spherical geometry. Geochem. Geophys. Geosystems 2004, 5. [Google Scholar] [CrossRef]
 Chan, F.K.; O'Neill, E.M. Feasibility study of a quadrilateralized spherical cube earth data base, Tech. Report EPRF 275 (CSC), 1975, Naval Environmental Prediction Research Facility.
 O’Neill, E.M.; Laubscher, R.E. Extended studies of a quadrilateralized spherical cube earth data base. Tech. Report NEPRF 376 (CSC), 1976, Naval Environmental Prediction Research Facility.
 Górski, K.M.; Hivon, E.; Banday, A.J.; Wandelt, B.D.; Hansen, F.K.; Reinecke, M.; Bartelmann, M. HEALPix: A framework for highresolution discretization and fast analysis of data distributed on the sphere. The Astrophysical Journal (ApJ) 2005, 622, 759–771. [Google Scholar] [CrossRef]
 Gibb, R.G. The rHEALPix Discrete Global Grid System. IOP Conf. Ser. Earth Environ. Sci. 2016, 34. [Google Scholar] [CrossRef]
 Bowater, D.; Stefanakis, E. 2018. The rHEALPix Discrete Global Grid System: Considerations for Canada. Geomatica 2018, 72, 27–37. [Google Scholar] [CrossRef]
 Béjar, R.; Lacasta, J.; LopezPellicer, F.J.; NoguerasIso, J. Discrete Global Grid Systems with quadrangular cells as reference frameworks for the current generation of Earth observation data cubes. Environmental Modelling & Software 2023, 162. [Google Scholar] [CrossRef]
 Lambers, M. Survey of Cube Mapping Methods in Interactive Computer Graphics. Vis. Comput. 2020, 36, 1043–1051. [Google Scholar] [CrossRef]
 Dimitrijević, A.; Lambers, M.; Rančić, D. Comparison of Spherical Cube Map Projections Used in PlanetSized Terrain Rendering. Facta Univ. Ser. Math. Inform. 2016, 31, 259–297. [Google Scholar]
 Lerbour, R.; Marvie, J.E.; Gautron, P. Adaptive realtime rendering of planetary terrains, in Full Paper Proc. Int. Conf. Computer Graphics, Visualization and Computer Vision (WSCG), 2010, 89–96.
 Grimm, C.; Niebruegge, B. Continuous Cube Mapping. J. Graph. GPU Game Tools 2007, 12, 25–34. [Google Scholar] [CrossRef]
 Nowell, P. Mapping a Cube to a Sphere. Available online: http://mathproofs.blogspot.de/2005/07/mappingcubetosphere.html (accessed on 24 May 2023).
 Purss, M. (Ed.) Topic 21: Discrete Global Grid Systems Abstract Specification, OGC 15104r5, 2017. Available online: http://docs.ogc.org/as/15104r5/15104r5.html (accessed on 24 May 2023).
 ISO 191701:2021, Geographic information — Discrete Global Grid Systems Specifications—Part 1: Core Reference System and Operations, and Equal Area Earth Reference System. 2021. Available online: https://www.iso.org/standard/32588.html (accessed on 24 May 2023).
 Gibb, R. (Ed.) Topic 21Discrete Global Grid SystemsPart 1 Core Reference system and Operations and Equal Area Earth Reference System”, OGC 20040r3. 2021. Available online: https://docs.ogc.org/as/20040r3/20040r3.html (accessed on 24 May 2023).
 Kmoch, A.; Matsibora, O.; Vasilyev, I.; Uuemaa, E. Applied opensource Discrete Global Grid Systems. AGILE GIScience Ser. 2022, 3, 41. [Google Scholar] [CrossRef]
 Purss, M.B.J.; Peterson, P.R.; Strobl, P.; Dow, C.; Sabeur, Z.A.; Gibb, R.G.; Ben, J. Datacubes: A Discrete Global Grid Systems Perspective. Cartogr. Int. J. Geogr. Inf. Geovisualization 2019, 54, 63–71. [Google Scholar] [CrossRef]
 Alderson, T.; Purss, M.; Du, X.; MahdaviAmiri, A.; Samavati, F. Digital Earth Platforms. In Manual of Digital Earth, 1st ed.; Guo, H., Goodchild, M.F., Annoni, A., Eds.; Springer: Singapore, 2020; pp. 25–54. [Google Scholar] [CrossRef]
 Digital Earth. Available online: https://jointresearchcentre.ec.europa.eu/scientificactivitiesz/digitalearth_en (accessed on 24 May 2023).
 Department of Defense, World Geodetic System 1984: Its Definition and Relationships with Local Geodetic Systems, Ver. 1.0.0, NGA.STND.0036_1.0.0_WGS84. 2014. Available online: https://nsgreg.nga.mil/doc/view?i=4085 (accessed on 24 May 2023).
 Petit, G.; Luzum, B. (Eds.), Terrestrial reference systems and frames, Technical Note No.36 in: IERS Conventions, 2010, 31–42. Available online: https://iersconventions.obspm.fr (accessed on 24 May 2023).
 Altamimi, Z.; Rebischung, P.; Métivier, L.; Métivier, L.; Collilieux, X. ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions. J. Geophys. Res. : Solid Earth 2016, 121, 6109–6131. [Google Scholar] [CrossRef]
 Adams, O. Latitude developments connected with geodesy and cartography: With tables including a table for Lambert equalarea meridional projection, U.S. Coast and Geodetic Survey Spec. Pub. No. 67, 1921.
 Snyder, J.P. Map Projections: A Working Manual; US Geological Survey Professional Paper No.1395; US Government Printing Office: Washington, DC, USA, 1987. [Google Scholar]
 Dimitrijević, A.; Strobl, P.; Lambers, M.; Milosavljević, A.; Rančić, D. Distortion Optimized Spherical Cube Mapping for Discrete Global Grid Systems. In ICIST 2020 Proceedings; Zdravković, M., Konjović, Z., Trajanović, M., Eds.; 2020; pp. 109–113. [Google Scholar]
 White, D.; Kimerling, A.J.; Overton, S.W. Cartographic and geometric components of a global sampling design for environmental monitoring. Cartogr. Geogr. Inf. Syst. 1992, 19, 5–22. [Google Scholar] [CrossRef]
 Greensky, J.B.S. G.; Czech, W.W.; Yuen, D.A.; Knox, M.R.; Damon, M.R.; Chen, S.S.; Kameyama, M.C. Ubiquitous interactive visualization of 3D mantle convection using a webportal with Java and Ajax framework. Vis. Geosci. 2008, 13, 105–115. [Google Scholar] [CrossRef]
 Chen, C.; Li, X.; Shen, X.; Xiao, F. Global shallow water models based on multimoment constrained finite volume method and three quasiuniform spherical grids. J. Comput. Phys. 2014, 271, 191–223. [Google Scholar] [CrossRef]
 Wongwathanarat, A.; Müller, E.; Janka, H.T. Threedimensional simulations of corecollapse supernovae: From shock revival to shock breakout. Astron. Astrophys. 2015, 577, A48. [Google Scholar] [CrossRef]
 Hara, K.; Inoue, K.; Urahama, K. Gradient operators for feature extraction from omnidirectional panoramic images. Pattern Recognit. Lett. 2015, 54, 89–96. [Google Scholar] [CrossRef]
 Barnes, R. Optimal orientations of discrete global grids and the Poles of Inaccessibility. Int. J. Digit. Earth 2019. [Google Scholar] [CrossRef]
 Reichard, C.G. Atlas des ganzen Erdkreises in der CentralProjection auf fechs Tafeln entworfen, Weimar. 1803. Available online: https://haabdigital.klassikstiftung.de/viewer/image/926414453 (accessed on 24 May 2023).
 Zhou, J.; Ben, J.; Wang, R.; Zheng, M.; Yao, X.; Du, L. A novel method of determining the optimal polyhedral orientation for discrete global grid systems applicable to regionalscale areas of interest. Int. J. Digit. Earth 2020, 13, 1553–1569. [Google Scholar] [CrossRef]
 Fuller, R.B. SynergeticsExplorations in the Geometry of Thinking; MacMillan Publishing Co. Inc.: New York, NY, USA, 1975. [Google Scholar]
 Sandvik, B. World Borders Dataset ver. 0.3. Available online: https://thematicmapping.org/downloads/world_borders.php (accessed on 24 May 2023).
Figure 1.
Unfolded regular polyhedra used in GDDGSs: (a) tetrahedron; (b) hexahedron; (c) octahedron; (d) dodecahedron; (e) icosahedron.
Figure 1.
Unfolded regular polyhedra used in GDDGSs: (a) tetrahedron; (b) hexahedron; (c) octahedron; (d) dodecahedron; (e) icosahedron.
Figure 2.
The process of transforming the planet surface into an addressable system of cells. The ellipsoid (a) is mapped onto the sphere (b), and the sphere onto the selected projection surfaces. The projection surfaces are usually the sides of circumscribed regular polyhedra or parts of cylinders (c). A partition is a projection of a part of the planet onto a projection surface. If the projection surfaces are not flat, they are unfolded (d) and finally divided into sections (e).
Figure 2.
The process of transforming the planet surface into an addressable system of cells. The ellipsoid (a) is mapped onto the sphere (b), and the sphere onto the selected projection surfaces. The projection surfaces are usually the sides of circumscribed regular polyhedra or parts of cylinders (c). A partition is a projection of a part of the planet onto a projection surface. If the projection surfaces are not flat, they are unfolded (d) and finally divided into sections (e).
Figure 3.
Deviations: (a) Auxiliary latitudes from geodetic latitude (θφ, θχ, and θβ); (b) Conformal from geocentric (χφ); (c) Approximated authalic from authalic (ββ’).
Figure 3.
Deviations: (a) Auxiliary latitudes from geodetic latitude (θφ, θχ, and θβ); (b) Conformal from geocentric (χφ); (c) Approximated authalic from authalic (ββ’).
Figure 4.
Two complementary partitions P0 (a) and P1 (b) of the DOEC projection, in the basic cylinder orientation and with no overlapping areas. The extent of the land mass and the graticule are displayed.
Figure 4.
Two complementary partitions P0 (a) and P1 (b) of the DOEC projection, in the basic cylinder orientation and with no overlapping areas. The extent of the land mass and the graticule are displayed.
Figure 5.
The distribution of angular, areal, and aspect distortion for both DOEC partitions, for the basic orientation of the projection cylinders.
Figure 5.
The distribution of angular, areal, and aspect distortion for both DOEC partitions, for the basic orientation of the projection cylinders.
Figure 6.
Rasterized world map without Antarctica in LatLon WGS 84 (EPSG:4326) projection used to test land mass distortion.
Figure 6.
Rasterized world map without Antarctica in LatLon WGS 84 (EPSG:4326) projection used to test land mass distortion.
Figure 7.
Two complementary partitions P0 (a) and P1 (b) of the DOEC projection obtained for optimally rotated projection cylinders (ϕ_{r} = 125°, θ_{r} = 50°, and ρ_{r} = 15°) to minimize landmass distortions and continental ruptures.
Figure 7.
Two complementary partitions P0 (a) and P1 (b) of the DOEC projection obtained for optimally rotated projection cylinders (ϕ_{r} = 125°, θ_{r} = 50°, and ρ_{r} = 15°) to minimize landmass distortions and continental ruptures.
Figure 8.
Graphical comparison of projections whose distortion parameters are listed in
Table 1:
(a) QSC;
(b) rHEALPix
^{E};
(c) rHEALPix
^{P};
(d) ASC;
(e) CCM;
(f) CSC;
(g) DOEC; and
(h) The threedimensional shape of the continents and the graticule.
Figure 8.
Graphical comparison of projections whose distortion parameters are listed in
Table 1:
(a) QSC;
(b) rHEALPix
^{E};
(c) rHEALPix
^{P};
(d) ASC;
(e) CCM;
(f) CSC;
(g) DOEC; and
(h) The threedimensional shape of the continents and the graticule.
Figure 9.
Comparison of the appearance of the partitions P0 for the cases of projection cylinders rotations: (a) R(125,50,15); (b) R(131,49,20).
Figure 9.
Comparison of the appearance of the partitions P0 for the cases of projection cylinders rotations: (a) R(125,50,15); (b) R(131,49,20).
Table 1.
The comparison of angular (ω), areal (σ) and aspect (α) distortions for the following projections: Dual Orthogonal Equidistant Cylindrical (DOEC), Adjusted Spherical Cube (ASC), Continuous Cube Mapping (CCM), Cartesian Spherical Cube (CSC), revised Hierarchical Equal Area isoLatitude Pixelization (rHEALPix), and Quadrilateralized Sphercal Cube (QSC). rHEALPix projection parameters are displayed independently for equatorial (rHEALPix^{E}) and polar (rHEALPix^{P}) regions. In addition to the minimum (min), maximum (max), and average (ave) values of the distortion parameters, the ratio between maximum and minimum areal distortion (σ_{max/min}), the ratio between average and minimum areal distortion (σ_{ave/min}), and the geometric mean of ω_{ave} and σ_{ave/min} (GM_{ωσ}) are displayed.
Table 1.
The comparison of angular (ω), areal (σ) and aspect (α) distortions for the following projections: Dual Orthogonal Equidistant Cylindrical (DOEC), Adjusted Spherical Cube (ASC), Continuous Cube Mapping (CCM), Cartesian Spherical Cube (CSC), revised Hierarchical Equal Area isoLatitude Pixelization (rHEALPix), and Quadrilateralized Sphercal Cube (QSC). rHEALPix projection parameters are displayed independently for equatorial (rHEALPix^{E}) and polar (rHEALPix^{P}) regions. In addition to the minimum (min), maximum (max), and average (ave) values of the distortion parameters, the ratio between maximum and minimum areal distortion (σ_{max/min}), the ratio between average and minimum areal distortion (σ_{ave/min}), and the geometric mean of ω_{ave} and σ_{ave/min} (GM_{ωσ}) are displayed.
Projection 
Angular distortion 
Areal distortion 
Aspect distortion 
GM_{ωσ}

ω_{min}

ω_{max}

ω_{ave}

σ_{min}

σ_{max}

σ_{ave}

σ_{max/min}

σ_{ave/min}

α_{min}

α_{max}

α_{ave}

QSC 
0.0 
25.081 
16.129 
1.910 
1.910 
1.910 
1.0 
1.0 
1.0 
1.555 
1.331 
4.016 
rHEALPix^{E}

0.0 
24.107 
7.964 
1.910 
1.910 
1.910 
1.0 
1.0 
1.0 
1.528 
1.155 
2.822 
rHEALPix^{P}

13.807 
49.250 
31.320 
1.910 
1.910 
1.910 
1.0 
1.0 
1.273 
2.429 
1.770 
5.596 
ASC 
0.0 
31.084 
11.572 
1.621 
2.293 
1.925 
1.414 
1.187 
1.0 
1.732 
1.234 
3.706 
CCM 
0.0 
31.084 
9.078 
1.463 
3.048 
1.967 
2.083 
1.344 
1.0 
1.732 
1.179 
3.493 
CSC 
0.0 
31.087 
11.489 
1.732 
2.309 
1.912 
1.333 
1.104 
1.0 
1.732 
1.235 
3.561 
DOEC 
0.0 
19.759 
5.864 
1.621 
2.293 
1.805 
1.414 
1.113 
1.0 
1.414 
1.113 
2.555 
Table 2.
The comparison of the averaged angular (ω_{ave}), areal (σ_{ave}), normalized areal (σ_{ave/min}) and aspect (α_{ave}) distortions and the geometric mean (GM_{ωσ}) for the initial orientation (R(0,0,0)), the proposed optimal rotation (R(125,50,15)), and the projection cylinder orientation causing the minimum landmass distortion (R(131,49,20)).
Table 2.
The comparison of the averaged angular (ω_{ave}), areal (σ_{ave}), normalized areal (σ_{ave/min}) and aspect (α_{ave}) distortions and the geometric mean (GM_{ωσ}) for the initial orientation (R(0,0,0)), the proposed optimal rotation (R(125,50,15)), and the projection cylinder orientation causing the minimum landmass distortion (R(131,49,20)).
R(ϕ_{r}, θ_{r}, ρ_{r}) 
ω_{ave}

σ_{ave}

σ_{ave/min}

α_{ave}

GM_{ωσ}

R(0,0,0) 
6.721 
1.832 
1.130 
1.130 
2.756 
R(125,50,15) 
3.557 
1.730 
1.067 
1.067 
1.948 
R(131,49,20) 
3.523 
1.729 
1.067 
1.066 
1.939 
Table 3.
Comparison of the effects of applying different auxiliary latitudes on the distortion of the P0 partition.
Table 3.
Comparison of the effects of applying different auxiliary latitudes on the distortion of the P0 partition.
Latitude 
Angular distortion 
Areal distortion 
Aspect distortion 
ω_{max}

ω_{ave}

σ_{min}

σ_{max}

σ_{max/min}

σ_{ave/min}

α_{max}

α_{ave}

Sphere 
19.758564 
5.866394 
1.621139 
2.292637 
1.414214 
1.113487 
1.414214 
1.113487 
Conformal 
19.758564 
5.866394 
1.621139 
2.300354 
1.418974 
1.114879 
1.414214 
1.113487 
Geocentric 
19.758882 
5.866548 
1.621139 
2.300349 
1.418971 
1.114877 
1.414222 
1.113490 
Authalic 
19.695547 
5.775854 
1.624769 
2.297771 
1.414214 
1.113487 
1.412636 
1.111732 
Approx. authalic 
19.695632 
5.775860 
1.624772 
2.297767 
1.414209 
1.113485 
1.412638 
1.111732 

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).