Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

134

Views

69

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

25 May 2023

Posted:

30 May 2023

You are already at the latest version

Alerts
Abstract
The rapid growth of Earth's global geospatial data necessitates an efficient system for organizing the data, facilitating data fusion from diverse sources, and promoting interoperability. Mapping the spheroidal surface of the planet presents significant challenges, as it involves balancing distortion and splitting the surface into multiple partitions. The distortion decreases as the number of partitions increases, but at the same time the complexity of data processing increases, since each partition represents a separate data set and is defined in its own local coordinate system. In this paper, we propose the dual orthogonal equidistant cylindrical projection method to mitigate distortion and reduce the number of partitions. Additionally, we use the rotation of projection cylinders to effectively minimize average angular and areal distortions of the Earth’s landmass and reduce the interruption of continental plates caused by partition edges. By incorporating auxiliary latitudes and proposing an approximate authalic latitude, we further enhance the mapping of the ellipsoid onto the sphere, simplifying calculations. Experimental results demonstrate a substantial reduction in distortion and interruption of continental plates. With only two partitions, an average landmass angular distortion of less than 3.56 degrees and an average normalized surface distortion of less than 1.07 were achieved.
Keywords: 
Subject: Computer Science and Mathematics  -   Information Systems

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 12N2 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.
φ = arctan ( ( 1 − e 2 ) ⋅ tan ( θ ) )
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).
χ = 2 ⋅ arctan ( tan ( π 4 + θ 2 ) ⋅ ( 1 − e ⋅ sin ( θ ) 1 + e ⋅ sin ( θ ) ) e 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 ci can be found in [29].
θ = χ + ∑ i = 1 ∞ c i ⋅ sin ( 2 χ )
The authalic latitude is calculated from the geodetic latitude using equations (4-6).
q = ( 1 − e 2 ) ⋅ ( sin ( θ ) 1 − e 2 ⋅ sin 2 ( θ ) − 1 2 ⋅ e ln ( 1 − e ⋅ sin ( θ ) 1 + e ⋅ sin ( θ ) ) )
q p = q θ = 90 ° = ( 1 − e 2 ) ⋅ ( 1 1 − e 2 − 1 2 ⋅ e ln ( 1 − e 1 + e ) )
β = arcsin ( q q p )
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:
β ' = arctan ( ( 1 − e 2 ) k ⋅ tan ( θ ) )
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 (180th 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).
θ P 1 = arcsin ( − cos ( θ P 0 ) ⋅ sin ( ϕ P 0 ) )
ϕ P 1 = − sgn ( θ P 0 ) ⋅ arccos ( − cos ( θ P 0 ) ⋅ sin ( ϕ P 0 ) / cos ( θ P 1 ) )
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 π 4 < ϕ p < 3 π 4 ) ∧ ( − π 4 < θ q < π 4 )
( ϕ p < − π 2 ∧ θ q < π 4 ) ∨ ( ϕ p > π 2 ∧ θ q > − π 4 )

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].
ω = 2 ⋅ arcsin ( | a − b | a + b )
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 rHEALPixE and rHEALPixP, respectively, in Table 1. Both the QSC and rHEALPixP 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 rHEALPixP 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 ω σ = ω a v e ⋅ σ a v e / min
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.

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 451-03-47/2023-01/200102.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Vestine, E.H.; Sibley, W.L.; Kern, J.W.; Carlstedt, J.L. Integral and Spherical-Harmonic Analyses of the Geomagnetic Field for 1955.0, Part 2. J. Geomagn. Geoelectr. 1963, 15, 73–89. [Google Scholar] [CrossRef]
  2. 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]
  3. 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]
  4. Sahr, K.; White, D.; Kimerling, A.J. Geodesic Discrete Global Grid Systems. Cartogr Geogr Inf Sci 2003, 30, 121–134. [Google Scholar] [CrossRef]
  5. Dimitrijević, A.; Rančić, D. Ellipsoidal Clipmaps–A planet-sized terrain rendering algorithm. Computers & Graphics 2015, 52, 43–61. [Google Scholar] [CrossRef]
  6. Kageyama, A.; Sato, T. Yin-Yang grid: An overset grid in spherical geometry. Geochem. Geophys. Geosystems 2004, 5. [Google Scholar] [CrossRef]
  7. Chan, F.K.; O'Neill, E.M. Feasibility study of a quadrilateralized spherical cube earth data base, Tech. Report EPRF 2-75 (CSC), 1975, Naval Environmental Prediction Research Facility.
  8. O’Neill, E.M.; Laubscher, R.E. Extended studies of a quadrilateralized spherical cube earth data base. Tech. Report NEPRF 3-76 (CSC), 1976, Naval Environmental Prediction Research Facility.
  9. Górski, K.M.; Hivon, E.; Banday, A.J.; Wandelt, B.D.; Hansen, F.K.; Reinecke, M.; Bartelmann, M. HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere. The Astrophysical Journal (ApJ) 2005, 622, 759–771. [Google Scholar] [CrossRef]
  10. Gibb, R.G. The rHEALPix Discrete Global Grid System. IOP Conf. Ser. Earth Environ. Sci. 2016, 34. [Google Scholar] [CrossRef]
  11. Bowater, D.; Stefanakis, E. 2018. The rHEALPix Discrete Global Grid System: Considerations for Canada. Geomatica 2018, 72, 27–37. [Google Scholar] [CrossRef]
  12. Béjar, R.; Lacasta, J.; Lopez-Pellicer, F.J.; Nogueras-Iso, 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]
  13. Lambers, M. Survey of Cube Mapping Methods in Interactive Computer Graphics. Vis. Comput. 2020, 36, 1043–1051. [Google Scholar] [CrossRef]
  14. Dimitrijević, A.; Lambers, M.; Rančić, D. Comparison of Spherical Cube Map Projections Used in Planet-Sized Terrain Rendering. Facta Univ. Ser. Math. Inform. 2016, 31, 259–297. [Google Scholar]
  15. Lerbour, R.; Marvie, J.-E.; Gautron, P. Adaptive real-time rendering of planetary terrains, in Full Paper Proc. Int. Conf. Computer Graphics, Visualization and Computer Vision (WSCG), 2010, 89–96.
  16. Grimm, C.; Niebruegge, B. Continuous Cube Mapping. J. Graph. GPU Game Tools 2007, 12, 25–34. [Google Scholar] [CrossRef]
  17. Nowell, P. Mapping a Cube to a Sphere. Available online: http://mathproofs.blogspot.de/2005/07/mapping-cube-to-sphere.html (accessed on 24 May 2023).
  18. Purss, M. (Ed.) Topic 21: Discrete Global Grid Systems Abstract Specification, OGC 15-104r5, 2017. Available online: http://docs.ogc.org/as/15-104r5/15-104r5.html (accessed on 24 May 2023).
  19. ISO 19170-1: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).
  20. Gibb, R. (Ed.) Topic 21-Discrete Global Grid Systems-Part 1 Core Reference system and Operations and Equal Area Earth Reference System”, OGC 20-040r3. 2021. Available online: https://docs.ogc.org/as/20-040r3/20-040r3.html (accessed on 24 May 2023).
  21. Kmoch, A.; Matsibora, O.; Vasilyev, I.; Uuemaa, E. Applied open-source Discrete Global Grid Systems. AGILE GIScience Ser. 2022, 3, 41. [Google Scholar] [CrossRef]
  22. 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]
  23. Alderson, T.; Purss, M.; Du, X.; Mahdavi-Amiri, 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]
  24. Digital Earth. Available online: https://joint-research-centre.ec.europa.eu/scientific-activities-z/digital-earth_en (accessed on 24 May 2023).
  25. 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).
  26. Petit, G.; Luzum, B. (Eds.), Terrestrial reference systems and frames, Technical Note No.36 in: IERS Conventions, 2010, 31–42. Available online: https://iers-conventions.obspm.fr (accessed on 24 May 2023).
  27. Altamimi, Z.; Rebischung, P.; Métivier, L.; Métivier, L.; Collilieux, X. ITRF2014: A new release of the International Terrestrial Reference Frame modeling non-linear station motions. J. Geophys. Res. : Solid Earth 2016, 121, 6109–6131. [Google Scholar] [CrossRef]
  28. Adams, O. Latitude developments connected with geodesy and cartography: With tables including a table for Lambert equal-area meridional projection, U.S. Coast and Geodetic Survey Spec. Pub. No. 67, 1921.
  29. 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]
  30. 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]
  31. 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]
  32. 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 web-portal with Java and Ajax framework. Vis. Geosci. 2008, 13, 105–115. [Google Scholar] [CrossRef]
  33. Chen, C.; Li, X.; Shen, X.; Xiao, F. Global shallow water models based on multi-moment constrained finite volume method and three quasi-uniform spherical grids. J. Comput. Phys. 2014, 271, 191–223. [Google Scholar] [CrossRef]
  34. Wongwathanarat, A.; Müller, E.; Janka, H.-T. Three-dimensional simulations of core-collapse supernovae: From shock revival to shock breakout. Astron. Astrophys. 2015, 577, A48. [Google Scholar] [CrossRef]
  35. 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]
  36. Barnes, R. Optimal orientations of discrete global grids and the Poles of Inaccessibility. Int. J. Digit. Earth 2019. [Google Scholar] [CrossRef]
  37. Reichard, C.G. Atlas des ganzen Erdkreises in der Central-Projection auf fechs Tafeln entworfen, Weimar. 1803. Available online: https://haab-digital.klassik-stiftung.de/viewer/image/926414453 (accessed on 24 May 2023).
  38. 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 regional-scale areas of interest. Int. J. Digit. Earth 2020, 13, 1553–1569. [Google Scholar] [CrossRef]
  39. Fuller, R.B. Synergetics-Explorations in the Geometry of Thinking; MacMillan Publishing Co. Inc.: New York, NY, USA, 1975. [Google Scholar]
  40. Sandvik, B. World Borders Dataset ver. 0.3. Available online: https://thematicmapping.org/downloads/world_borders.php (accessed on 24 May 2023).
Preprints 74656 g001
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.
Preprints 74656 g001
Preprints 74656 g002
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).
Preprints 74656 g002
Preprints 74656 g003
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 (β-β’).
Preprints 74656 g003
Preprints 74656 g004
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.
Preprints 74656 g004
Preprints 74656 g005
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.
Preprints 74656 g005
Preprints 74656 g006
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.
Preprints 74656 g006
Preprints 74656 g007
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.
Preprints 74656 g007
Preprints 74656 g008
Figure 8. Graphical comparison of projections whose distortion parameters are listed in Table 1: (a) QSC; (b) rHEALPixE; (c) rHEALPixP; (d) ASC; (e) CCM; (f) CSC; (g) DOEC; and (h) The three-dimensional shape of the continents and the graticule.
Figure 8. Graphical comparison of projections whose distortion parameters are listed in Table 1: (a) QSC; (b) rHEALPixE; (c) rHEALPixP; (d) ASC; (e) CCM; (f) CSC; (g) DOEC; and (h) The three-dimensional shape of the continents and the graticule.
Preprints 74656 g008
Preprints 74656 g009
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).
Preprints 74656 g009
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 (rHEALPixE) and polar (rHEALPixP) 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 (rHEALPixE) and polar (rHEALPixP) 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
rHEALPixE 0.0 24.107 7.964 1.910 1.910 1.910 1.0 1.0 1.0 1.528 1.155 2.822
rHEALPixP 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.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated