Efficient Distortion Mitigation and Partition Reduction in Mapping Earth's Geospatial Data: Dual Orthogonal Equidistant Cylindrical Projection Approach

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 high-resolution 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 re-projections 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 mid-20th 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 semi-regular circumscribed polyhedra. Commonly used base polyhedra are the five Platonic solids, in particular the hexahedron (cube), while the most common semi-regular 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 out-of-core planet-scale 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 Yin-Yang grid [6], as projection planes to reduce the partitions to two. While the Yin-Yang 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² points as uniformly as possible over the surface of the unit sphere [9]. These hybrid projections combine the Lambert cylindrical equal-area 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 equal-area, 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 planet-scale 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 19170-1: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 Yin-Yang 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 two-dimensional 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 equal-area projections preserve area but distort shape considerably. A projection that is both conformal and equal-area 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 Geospatial-Intelligence 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.

$$\phi =\arctan \left((1-e^2)\cdot \tan (\theta )\right)$$

(1)

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 equal-area mapping, eliminating areal distortion. The conformal latitude can be computed from the geodetic latitude using equation (2).

$$\chi =2\cdot \arctan \left(\tan \left(\frac{\pi }{4}+\frac{\theta }{2}\right)\cdot {\left(\frac{1-e\cdot \sin (\theta )}{1+e\cdot \sin (\theta )}\right)}^{\frac{e}{2}}\right)-\frac{\pi }{2}$$

(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.0E-12. The corresponding values for the coefficients c_i can be found in [29].

$$\theta =\chi +{\displaystyle \sum _{i=1}^{\infty }{c}_i\cdot \sin (2\chi )}$$

(3)

The authalic latitude is calculated from the geodetic latitude using equations (4-6).

$$q=\left(1-e^2\right)\cdot \left(\frac{\sin (\theta )}{1-e^2\cdot {\sin }^2(\theta )}-\frac{1}{2\cdot e}\ln \left(\frac{1-e\cdot \sin (\theta )}{1+e\cdot \sin (\theta )}\right)\right)$$

(4)

$$q_p=q_{\theta =90°}=\left(1-e^2\right)\cdot \left(\frac{1}{1-e^2}-\frac{1}{2\cdot e}\ln \left(\frac{1-e}{1+e}\right)\right)$$

(5)

$$\beta =\arcsin \left(\frac{q}{q_p}\right)$$

(6)

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.4E-2° (50.4"). Considering the small deviation from conformal latitude, the relative simplicity of the calculation, and the availability of a straightforward closed-form 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 closed-form inverse transformation, we propose the following approximated formula:

$$\beta \text{'}=\arctan \left({(1-e^2)}^k\cdot \tan (\theta )\right)$$

(7)

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.16E-5° (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 semi-regular 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 semi-regular 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 Yin-Yang [6] is obtained. The Yin-Yang 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 core-collapse 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 Yin-Yang 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 anti-meridian (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).

$${\theta }_{P1}=\arcsin \left(-\cos ({\theta }_{P0})\cdot \sin ({\varphi }_{P0})\right)$$

(8)

$${\varphi }_{P1}=-\mathrm{sgn}({\theta }_{P0})\cdot \arccos \left(-\cos ({\theta }_{P0})\cdot \sin ({\varphi }_{P0})/\cos ({\theta }_{P1})\right)$$

(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).

$$\left(-\frac{3\pi }{4}<{\varphi }_p<\frac{3\pi }{4}\right)\wedge \left(-\frac{\pi }{4}<{\theta }_q<\frac{\pi }{4}\right)$$

(10)

$$\left({\varphi }_p<-\frac{\pi }{2}\wedge {\theta }_q<\frac{\pi }{4}\right)\vee \left({\varphi }_p>\frac{\pi }{2}\wedge {\theta }_q>-\frac{\pi }{4}\right)$$

(11)

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) semi-axes of the indicatrix. The formulas for calculating the indicatrices and the corresponding deformation parameters can be found in [29].

$$\omega =2\cdot \arcsin \left(\frac{\left|a-b\right|}{a+b}\right)$$

(12)

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 semi-axes 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 in-depth 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 location-specific 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 on-the-fly 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 out-of-core 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 high-resolution aerial imagery at multiple levels of detail. For out-of-core 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 three-dimensional 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 equal-area 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 equal-area 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).

$$G{M}_{\omega \sigma }=\sqrt{{\omega }_{ave}\cdot {\sigma }_{ave/\min }}$$

(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 world-scale 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.