Regression and Evaluation on a Forward Interpolated Version of the Great Circle Arcs–Based Distortion Metric of Map Projections

: We studied the numerical approximation problem of distortion in map projections. Most widely used differential methods calculate area distortion and maximum angular distortion using partial derivatives of forward equations of map projections. However, in certain map projections, partial derivatives are difﬁcult to calculate because of the complicated forms of forward equations, e.g., equations with iterations, integrations, or multi-way branches. As an alternative, the spherical great circle arcs–based metric employs the inverse equations of map projections to transform sample points from the projection plane to the spherical surface, and then calculates a differential-independent distortion metric for the map projections. We introduce a novel forward interpolated version of the previous spherical great circle arcs–based metric, solely dependent on the forward equations of map projections. In our proposed numerical solution, a rational function–based regression is also devised and applied to our metric to obtain an approximate metric of angular distortion. The statistical and graphical results indicate that the errors of the proposed metric are fairly low, and a good numerical estimation with high correlation to the differential-based metric can be achieved.


Introduction
A geographic map is a planar symbolic depiction describing natural or social features and characteristics on Earth or other planets [1]. In the process of map making, map projections [2][3][4][5] are fundamental because they represent an accurate transformation between three-dimensional and planar coordinates.
As spheres and ellipsoids are not developable surfaces, distortions in map projections are inevitable [22,23]. The measurement and analyses of distortions [24,25] for different types of map projections are critical and still a challenging task for the design [16], categorization, evaluation, comparison, selection [26,27], and optimization [19] of map projections.
Nicolas Tissot devised a visual indicator (i.e., Tissot's indicatrix or Tissot's ellipse) to represent linear, angular, and areal distortion in map projections [28,29]. Many authors derived different computational formulae of Tissot's original theory [28] in the numerical details [24,30]. Among the different forms of equations, area distortion and maximum angular distortion are calculated by the partial derivatives of forward equations (see p. 24 in [2]), and they are used as the reference of our metric.
The singular value decomposition (SVD)-based method [30] employs an algebraic approach. The SVD-based method can alternatively be used to calculate the same results as the maximal and minimal linear scales at a given point. There are also several metrics proposed based on partial derivatives. For example, the large-scale distortions of flexion and skewness are defined as Goldberg-Gott indicatrices [31,32].
However, above the Airy-Kavrayskiy metric, Goldberg-Gott indicatrices, Capek Q, and the SVD-based methods are all dependent on partial derivatives of forward or inverse equations of map projections.
Partial derivatives of forward equations are not always yielded in convenient ways. For example, Mollweide, Boggs Eumorphic, Eckert IV, Nell, and Putnin . š P 2 projections do not have the closed-form of forward equations, and they require iterative computations [2,34], such as Newton-Raphson iterations. In addition, certain projections may be interrupted (e.g., Hierarchical Equal Area isoLatitude Pixelisation, or HEALPix [35,36], interrupted Goode homolosine [22,23,37] and icosahedral Snyder equal area [38] projections), derivative discontinuous (e.g., uninterrupted Goode homolosine [14] and Robinson [39] projections), have singularity points (Mercator and Peirce Quincuncial [40] projections), or may have an elliptic integral (Adams World in a Square I and II and Peirce Quincuncial [40] projections) in forward equations. Another complicated example is the Chamberlin trimetric projection [41], which has multiple loops and several if-else branches in forward equations in the implementation of PROJ library. The above projections and similar cases indicate that calculating partial derivatives, as well as the differential metric and its derivatives, for map projections is not always an easy task.
In addition to map projections defined using mathematical formulae, tessellation using direct spherical subdivision (DSS) produces discrete global grids [42,43]. DSS methods are recursive subdivisions; hence, there is no explicit form of transformations between the geographic coordinate system and grid coordinate system. Consequently, the differential methods and its derivatives are not applicable to DSS-based grids [44].
The averaged ratio metric between complementary profiles [44] is based on spherical small circle arcs (SCA). SCA-based metrics and their derived spherical great circle arc (GCA)-based indicators [45] are all differentially independent. These methods sample points on the projection plane and transform sample points from the projection plane to the sphere to calculate SCA-or GCA-based metrics. The primary benefit of SCA-or GCA-based metrics includes independence from differential-based calculation. However, one shortcoming is their dependency on inverse equations of map projections. Nonetheless, certain map projections do not have explicit forms of inverse equations, such as Apian Globular I, Bacon Globular, Ortelius oval, Denoyer semi-elliptical, Ginsburg VIII, Laskowski, and Larrivée projections. Implicit forms of inverse equations for these map projections are also difficult to derive. This is a restriction of SCA-or GCA-based metrics when evaluating map projections without explicit or implicit forms of inverse equations.
In this paper, we extend the inverse equation-dependent GCA-based metric [45] to a forward interpolated form that is solely dependent on the forward equations of map projections. The proposed metric is abbreviated as forward GCA (FWD-GCA) for conciseness. The FWD-GCA approach samples points on the sphere, and transforms these sampled points to a projection plane. Then, the vector interpolation is applied to the sample points on the sphere to obtain a new set of sample points to create the new projected points on the projection plane that satisfies the orthogonal conditions of the GCA-based metric.
Regression between the FWD-GCA metric and differential-based maximum angular distortion was also studied. For achieving a high precision of regression, rational functions with degree one were employed, and the relationship among regression coefficients was studied and determined.
The FWD-GCA metric is derived and implemented for the unit sphere in this paper. Under the assumption of the spherical Earth, one known shortcoming of the FWD-GCA metric is that it is currently unsuitable for ellipsoidal map projections.
The remainder of this paper is organized as follows. Relevant studies related to differential methods and GCA-based metrics are outlined in Section 2. The FWD-GCA metric, and the process of calculation and interpolation are introduced in Section 3. Regressions by employing linear and rational functions are established and presented in Section 4. The experiments and validation are shown in Section 5. In Section 6, selected map projections are evaluated and compared. Finally, Section 7 discusses the results and concludes this paper.

Differential Metric
Given partially differentiated forward equations x(ϕ, λ) and y(ϕ, λ) of a map projection, where ϕ and λ represent the latitude and longitude, respectively, and x and y represent the projected coordinates, the differential metric (see p. 24 in [2]) calculates the maximum angular distortion ω (see Equation (6)) and area distortion s (see Equation (9)) for the projections on the sphere. Metric ω represents the greatest deviation from the correct angle in Tissot's indicatrix. The implementation of the differential metric in Julia is provided in Codes S1 (for map projections with the explicit form of forward equations) and S2 (for map projections with the implicit form of forward equations) in supplementary materials. For the simplicity of presentation, the formulae for the differential metric (see Equations (1)-(9)), GCA-based metric (see Section 2.2), and FWD-GCA metric (see Section 3) are defined on a unit sphere S 2 centered at origin O.
where h and k represent scale factors along the meridian and the parallel, respectively; a and b represent the maximum and minimum linear scale factors, respectively; ω denotes the maximum angular distortion; and s represents the areal scale factor [2,25,32]. The maximum angular distortion ω and area distortion s are used as reference metrics (see Sections 4 and 5) for evaluating the GCA-based and FWD-GCA metrics in this paper.
Using differential-based metrics is straightforward for map projection with explicit forward equation (see Code S1). However, the Newton-Raphson iteration is needed for map projections with implicit forward equations (see Code S2); the derivation of partial derivatives is also required (see Appendix A).

Gca-Based Metric
The averaged ratio metric between complementary profiles was introduced in [44], and uses small circle arcs on a sphere to evaluate the shape regularity of tessellation grids or map projections. Instead of small circle arcs, the simplified GCA-based metric [45] uses great circle arcs (see black arcs in Figure 1b) to calculate the area distortion s gc (see Equation (22)) and angular distortion ρ gc (see Equation (20)). Ref. [45] revealed that the metric ρ gc is a linear estimation of the differential-based metric ρ d (see Equation (21)).
The GCA-based metric can be calculated if the inverse equations of map projection are available. Given a map projection, the forward equations are assumed to be x(ϕ, λ) and y(ϕ, λ), and the inverse equations are assumed to be ϕ(x, y) and λ(x, y), where ϕ and λ are latitude and longitude coordinates, respectively, and x and y are planar coordinates on the given two-dimensional projection plane. Consider a point P 0 = (X 0 , Y 0 , Z 0 ) T on the unit sphere S 2 ⊂ R 3 , the geographic coordinates of point P 0 can be calculated according to: where ϕ 0 and λ 0 represent the latitude and longitude coordinates of point P 0 , respectively. The GCA-based metric can be calculated as:  First, calculate the coordinates of the projected point p 0 (see Figure 1a) on the projection plane by using the forward equations of map projection.
Next, calculate the four corner points (see Figure 1a) of an infinitesimal square on the projection plane as: where i = 1, 2, 3, 4. We assume here that p 1 and p 3 are located on one diagonal of the square, p 2 and p 4 are located on the other diagonal of the square, and is a sufficiently small value ( = 10 −6 in this paper) representing half the length of the side of square. Points p i , i = 1, 2, 3, 4, are located on the corners of a square, which is an undistorted quadrilateral on a plane. The diagonals of the square act as an orthogonal frame on the projection plane for the GCA-based metric. Inspired by the GCA-based method, the proposed FWD-GCA method (see Section 3) uses an orthogonal frame on the projection plane to establish metrics for evaluating the distortions of the map projections. Now, let us consider the GCA-based method. We calculate the geographic coordinates (ϕ i , λ i ) T for each two-dimensional point p i using the inverse equations of the map projection as: where i = 1, 2, 3, 4. Subsequently, we calculate the Cartesian coordinates (i.e., the earth-centered, earthfixed (ECEF) Cartesian coordinates) of the three-dimensional points P i (see Figure 1b) on the unit sphere S 2 for each two-dimensional point p i , i = 1, 2, 3, 4.
where i = 1, 2, 3, 4. Now, prepare intermediate three-dimensional vectors q i (i = 1, 2, see Figure 1b) and normalized three-dimensional vectorsq i (i = 1, 2) for the GCA-based metric, and calculate r i (i = 1, 2, see Figure 1b), which is the length of the great circle arc connecting points P i and P i+2 on the sphere, where i = 1, 2.
where i (or i for vectors q i andq i ) = 1, 2. Next, calculate the length and angle components of the GCA-based metric.
Finally, combine the above components as a GCA-based metric by the arithmetic mean.
The GCA-based metric ρ gc has a strong (>0.988) linear correlation [44] to a differentialbased shape regularity metric ρ d .
where the maximum and minimum linear scale factors a and b are calculated according to Equations (7) and (8).
The reciprocal of ρ d was at first defined as the metric of angular distortion, which is equal to b/a [25,32,46]. The metric ρ d is then adopted in [44,45] to measure the shape distortion of the tessellation.
The area distortion is discussed in [45]. It is defined as the ratio between the area of the infinitesimal region on the projection plane and the area of the corresponding region on the sphere; s gc is introduced as an alternative to the differential-based areal scale factor s in Equation (9).
where 4 2 is the area of infinitesimal square on the projection plane (see Figure 1a).

Forward Interpolated Version of Gca-Based Metric
The calculation of the GCA-based metric is dependent on ϕ(x, y) and λ(x, y), the inverse equations of the map projection (see Equation (13)). To avoid the influence of inverse equations, the FWD-GCA metric, an interpolated version of the GCA-based metric, is introduced in this section. For a clear and intuitive comparison, Figure 2 summarize the principle of and difference among different methods. The GCA-based method generates four points, around the given point on the projection plane (R 2 ), while the proposed FWD-GCA-based method directly generates four points on the unit sphere S 2 ⊂ R 3 . The FWD-GCA-based metric has two main steps. In Section 3.1, the sample points on the unit sphere S 2 (on red lines in Figure 3a) are generated and converted to the projection plane (on red lines in Figure 3b). In Section 3.2, vector interpolation is applied to sample points on the projection plane to obtain orthogonal vectors (blue lines in Figure 3b), and the points on the orthogonal vectors (blue lines in Figure 3b) on the projection plane are converted to a sphere (denoted asp i on the blue lines in Figure 3c, i = 1, 2, 3, 4) to calculate the FWD-GCA metric. An implementation of the FWD-GCA metric on the unit sphere S 2 in Julia is provided in Code S3 in the supplementary materials.

Spherical Sampling and Forward Projection
Given a three-dimensional point P 0 = (X 0 , Y 0 , Z 0 ) T on the unit sphere S 2 , calculate the geographic coordinates (ϕ 0 , λ 0 ) T of point P 0 according to Equation (10) and the twodimensional projection coordinates p 0 = (x 0 , y 0 ) T according to Equation (11).
The calculation for the FWD-GCA metric is as described below. First, let E be a three-dimensional point on S 2 and near P 0 . For the ease of the calculation, E is offset along the z-axis (i.e., on the line between the North and South Poles, with positive values increasing northward) of ECEF Cartesian coordinates, whereas the singularity case is avoided (if Z 0 = ±1) using a sign function of Z 0 (see Equation (24)).
where is a sufficiently small value, and = 10 −6 in this paper. Next, calculate four unit three-dimensional vectors as e 1 = normalize(E − P 0 ) and e i+1 = p 0 × e i , where i = 1, 2, 3; the three-dimensional vector p 0 = P 0 − O is the vector to P 0 in R 3 ; and O denotes the origin of S 2 . Now, calculate the four points (see black dots in Figure 3a) along the orthogonal axes (see red lines in Figure 3a) on S 2 according to: where i (or i for vectors e i ) = 1, 2, 3, 4, and = 10 −6 in this paper. Convert the Cartesian coordinates of the sphere (Figure 3a) to the geographic coordinates Φ i and Λ i using Equation (26); calculate four basis points b i on the projection plane (see Figure 3b) from Φ i and Λ i , i = 1, 2, 3, 4, using the forward equations of map projection (see Equation (27). The uppercase Greek letters Φ and Λ with subscripts in Equations (26) and (27) are used to differentiate them with ϕ and λ with subscripts in Equations (13) and (14). Figure 4 demonstrates more details and annotation than Figure 3b. In Figure 4, there are different cases of the relative positions of the points b i depending on the specific map projection and the positions of given geographic coordinates ϕ 0 and λ 0 .
where i = 1, 2, 3, 4. On the projection plane (see Figure 3b or Figure 4), calculate four two-dimensional basis vectors as: where i (or i for vectors b i ) = 1, 2, 3, 4. In Figure 4, vectors b i and b i+2 appear to be collinear, where i = 1, 2. However, this pair of vectors is not strictly collinear, as in Equation (25) is not infinitely small.
Calculate the two normalized two-dimensional basis vectorsb i (marked as blue lines in Figure 4) as:b where i = 1, 2. Additionally, calculate the angle α shown in Figure 4.

Interpolation of Vectors on Sphere
Similar to sample points (see Figure 1a) on the projection plane coinciding with corners of a square in the GCA-based metric, the FWD-GCA-based metric needs to establish sample points on the projection plane to satisfy the relationship of orthogonal diagonals (see blue lines in Figure 3b). In Figure 4, red and yellow lines represent the candidate orthogonal frames for the FWD-GCA-based metric.
Interpolation based on the orthogonal frames is as described below. Firstly, establish and solve a set of equations (see Equations (32)- (35) or (40)-(43)) on the projection plane (see Figure 3b) to obtain the coefficients of interpolation depending on the sign of α. Then, the coefficients are used to interpolate the point (see Equations (36)- (39) or (44)-(47)) on the sphere to obtain P i , i = 1, 2, 3, 4 (see Figure 3c).
If α ≥ 0 (see Figures 4a), select b 1 and b 2 as the basis vector to solve Equations (32) and (33) and obtain interpolation coefficients for the two vectors (red lines in Figure 4a), and select b 3 and b 4 as the basis vector to solve Equations (34) and (35) and obtain interpolation coefficients for the other two vectors (yellow lines in Figure 4a).
The coefficient of interpolation s ij is calculated, where i = 1, 2, 3, 4 and j = 1, 2, 3, 4; we calculate points on sphere by vector interpolation as: If α < 0 (Figure 4b), select b 1 and b 4 as the basis vectors to solve Equations (40) and (43), respectively, and obtain interpolation coefficients for the two vectors (red lines in Figure 4b). Similarly, select b 2 and b 3 as the basis vector to solve Equations (41) and (42) and obtain interpolation coefficients for the other two vectors (yellow lines in Figure 4b).
The coefficient of interpolation s ij is calculated, where i = 1, 2, 3, 4 and j = 1, 2, 3, 4; we calculate points on the sphere by vector interpolation as: By substituting three-dimensional points P i (see Equations (36)- (39) or (44)-(47), and Figure 3c), i = 1, 2, 3, 4, into the right-hand side of Equation (15), we obtain threedimensional vectors q i , i = 1, 2. Next, we obtain a forward version of the GCA-based metric (or FWD-GCA metric) according to Equations (16) to (20); it is denoted as ρ f wd . Area distortion in the FWD-GCA metric can be obtained by rewriting Equation (22) as: where is half the length of the diagonal of an infinitesimal square on the projection plane ( is half the length of side of an infinitesimal square in Equation (22)), and 2 2 is the area of an infinitesimal square on the projection plane in Equation (48), while 4 2 is the area of an infinitesimal square on the projection plane in Equation (22).

Linear Regression
According to the results in [44,45], the differential-based ρ d and GCA-based metric ρ gc are highly correlated (> 0.988 for 14 non-conformal map projections [45]); hence, ρ gc can be approximately estimated from ρ d , and vice versa using the linear regression model of ordinary least squares (OLS). That is, where ρ gc d denotes a metric derived from metric ρ gc ; it is an approximation of metric ρ d . An example of metric distributions of ρ d and ρ gc is shown in Figure 5a,b with respective scalar bars for the Bonne projection with first standard parallel ϕ 1 = 30 • . The Pearson product-moment correlation coefficient (PPMCC) between ρ d and ρ gc for the Bonne projection is 0.9883, which is almost the worst result among the 14 non-conformal map projections in [45].
The black contours in Figure 5a are derived from the differential metric ρ d , and the red contours in Figure 5b are derived from the metric ρ (c) FWD-GCA metric ρ f wd (distribution) andρ f wd (contours). In Section 3, the proposed FWD-GCA metric is constructed from the GCA-based metric using the same mechanism of the GCA-based metric. As in the case of the GCAbased metric, we assume that the proposed FWD-GCA metric should also be a linear approximation of the differential metric ρ d , i.e., where ρ f wd d denotes a metric derived from the metric ρ f wd , and it should be an approximation to the metric ρ d .
OLS-based regression can be applied to Equation (50)

Rational Function Regression
Now, let us consider the maximum angular distortion ω in Equation (6). Using Equations (6)-(8) and (21), we derive that:  (49) and (50)) of ρ d by OLS-based regression, we substitute the right-hand sides of Equations (49) and (50) into Equation (51); we can further assume that: where coefficients β  As shown in Equations (52) and (53), ω is represented by different forms of rational functions of GCA-based metric ρ gc or FWD-GCA-based metric ρ f wd . Without the loss of generality, we rewrite Equations (52) and (53) and establish a general form for rational functions with degree one (or linear rational functions) for independent variables ρ gc or ρ f wd using unknown coefficients δ gc i or δ f wd i , where i = 0, 1, 2, as:

Approach of Estimating Regression Coefficients
Equation (55) is a proposed equation to predict ω from the FWD-GCA-based metric with unknown coefficients. Therefore, we need to further establish an approach to estimate the value of coefficients, which can covers all different types of map projections, e.g., conformal, equal-area or compromise map projections.
Firstly, we would compare linear regression with the proposed rational function-based regression (see Section 5.2.1) to find out the distribution and characteristic of regression coefficients by considering as many map projections as possible. However, there are hundreds of map projections and their variants; therefore, one cannot enumerate as many of them as possible. That means it is impossible to make a regression for all innumerable existing map projections.
A practical approach is to use a limited number of randomly selected map projections to obtain individual regression coefficients δ f wd i , i = 0, 1, 2, for each map projection. Next, these individual regression coefficients are used to obtain synthesized coefficients based on the distribution and characteristic of regression coefficients. As the synthesized coefficients are obtained and validated, we can use Equation (55) to calculate ω f wd .
To further eliminate the effect of individual map projection, we can use certain groups of randomly selected map projections (see Section 5.2.2) to obtain individual regression coefficients, for each group of map projections. Then, we can obtain synthesized coefficients from these individual regression coefficients.
As only a limited numbers of map projections are used for regression, it is necessary to study the sensitivity (see Section 5.3.3) of regression coefficients δ f wd i , i = 0, 1, 2, yielding different sources of the selected map projections.
Next, we will estimate the coefficients of the above rational functions with degree one in Equation (55), as well as validate the accuracy of the above rational function-based regression method in Section 5.

Overview of Implementation and Experiments
More than 80 map projections for small-scale world maps in PROJ library [20] were evaluated in our study on a laptop with 64-bit quad-core of 1.9 GHz CPU, 16 GB memory, and 64-bit Manjaro Linux 21.0.3. Julia [47] is a modern, high-level, high-performance, and dynamic programming language. Differential metric (see Codes S1 and S2 in supplementary materials, and see Appendix A), GCA-based metric, and FWD-GCA metric (see Code S3) were implemented in Julia because of its speed and flexibility [47]. However, certain calculations, analyses, and plots were implemented in Python owing to its simplicity and usability. The versions of Julia and Python used were 1.6.1 and 3.9.4, respectively. Two additional Julia packages, Proj4.jl and ForwardDiff.jl [48], were employed in our implementation. Proj4.jl is a wrapper of PROJ library, which is a coordinate conversion tool. It is used to transform geographic coordinates to projection coordinates, and vice versa. The underlying version of PROJ used in Proj4.jl is 7.2.1. ForwardDiff.jl is used to calculate partial derivatives, which is used in differential methods for comparisons. A few additional Python packages, such as numpy, scipy, healpy, matplotlib, were also used in our experiments.
For all map projections we evaluated, we calculated the differential-based metrics s and ρ d , and obtained ω by using Equation (51). We also calculated the GCA-based metrics s gc and ρ gc , and the proposed FWD-GCA metrics s f wd and ρ f wd . Then, ω gc and ω f wd were calculated for validating the approach of the proposed rational function-based regression. Approximately 200k sample points on the sphere were generated using healpy, a wrapper of HEALPix [35] algorithm, to evaluate the metrics s gc , ρ gc , ω gc , s f wd , ρ f wd , and ω f wd . Parameter N side = 128 was used in the HEALPix algorithm to generate approximately 200k (12 × N 2 side ) near-uniformly distributed sample points on the sphere for all map projections to be evaluated. All 200k sample points were shifted by 0.003 radians along the longitudinal direction to avoid λ = ±π.

Distribution of Regression Coefficients
Firstly, we compared the linear regression and rational function-based regression between the differential metric and FWD-GCA-based metric for more than 80 map projections in the PROJ library.      Figure 8a Most † Non-Cylindrical, Non-Conic, and Non-Azimuthal Projections in Figure 8b range of δ 0.999999999999992 0.9998 † Rectangular polyconic (see letter "A" in Figure 8a) is a non-cylindrical, non-conic, and non-azimuthal projection, but it is not in Figure 8b.

All Non-Conformal Map Projections in
Letters P to W in Figure 8b represent the same map projections as those marked as Letters A, B, and E-K in Figure 7 for comparison. It can be observed from Figures 7 and 8b that regression coefficients of these map projections deviate from the blue line in Figure 7 but coincide to blue or red lines in Figure 8b. That means that rational function-based regression can produce coefficients satisfying a highly correlated linear relationship (see PPMCC in Table 1).
Although regression coefficients δ f wd i in Figure 8a,b have a fairly large range, i = 0, 1, 2, we assumed that they should have an approximately linear relationship (see Equations (56) and (57)) with additional regression coefficients γ i , i = 0, 1, 2, 3, according to the results shown in Figure 8.
As linear relationships (see Equations (56) and (57)) are intuitively established for the rational function-based regression coefficients, we needed to further study the effectiveness of the above linear relationships as well to obtain definite regression coefficients δ f wd i (see Equation (55)), i = 0, 1, 2, to estimate the maximum angular distortion ω using Equation (55) from the proposed the FWD-GCA metric ρ f wd .

Regression Results of Random Selection of Map Projections
First, we studied a limited number of map projections whose rational function-based regression coefficient δ f wd 0 is less than 10 (see Figure 8b). We randomly selected 12 groups of 10 map projections from Figure Table S2. Table 2 summarizes Table S2. Table 2 indicates that the ranges of β

Candidate of Regression Coefficients
As rational function-based regression (see Equation (55)) is established, we need definite coefficients (i.e., δ f wd i in Equation (55), i = 0, 1, 2) for the regression. Now, let us select candidates for regression coefficients and study the sensitivity of the regression coefficients.
According to Table 2 Table 3. We also chose three pairs of β f wd 0 and β f wd 1 (see Table 3, two of them derived from Table S2 and one manually) for comparisons of linear regression. Rational function-based regression also uses a manual setting of coefficients.  Table 3 reveals that rational function-based regression produces the same or similar coefficients δ f wd i , i = 0, 1, 2, which are derived from γ f wd i (see Table S2), i = 0, 1, 2, 3. Any values the same or similar of δ f wd i , i = 0, 1, 2 in Table 3, were not selected for further evaluation. Regression coefficients used for further sensitivity analyses of linear regression and rational function-based regression are marked in the fourth and last columns of Table 3, respectively.

Sensitivity of Regression Coefficients for Linear Regression
We used candidates of linear regression coefficients in Table 3 to calculate ω f wd by using Equations (53) and (55). The distributions of ω f wd calculated by different linear regression coefficients, and the corresponding contours of ω f wd (in red) and ω (in black), are plotted in Fiugre 9 for Putnin . š P 4 , Winkel Tripel, Eckert II, and Collignon projections. For each map projection shown in Figure 9, we calculate ω f wd and ω for approximately one million sample points on the projection plane. Contours of the metrics were also generated based on one million sample points.
Due to the symmetry of Putnin . š P 4 , Winkel Tripel, and Eckert II projections, we only plotted a quarter area (0 ≤ ϕ ≤ π/2 and 0 ≤ λ ≤ π) on the projection plane in Figure 9a-c. The hemisphere of Collignon projection is plotted in Figure 9d due to the same reason. PPMCC and root-mean-square error (RMSE) (see Equations (58) and (59)) were used to evaluate the correlation and error between the proposed FWD-GCA metric-based ω f wd . We calculated PPMCC and RMSE (see values in  by using approximately 200k sample points on the sphere, which were generated using the HEALPix algorithm for the sake of consistency of evaluating different map projections.
PPMCC ω, where n is the number of sample points. According to the results in Figure 9, linear regression does not achieve a good estimation for all selections of β

Sensitivity of Regression Coefficients for Rational Function-Based Regression
Now, we study the sensitivity of regression coefficients δ f wd i , i = 0, 1, 2, for rational function-based regression. Distribution and contours based on the results of rational function-based regression are plotted using approximately one million sample points for Putnin . š P 4 (Figure 10), Winkel Tripel (Figure 11), Eckert II (Figure 12), and Collignon ( Figure 13)   According to Figures 10-13, using the proposed FWD-GCA metric with rational function-based regression to estimate the maximum angular distortion ω is insensitive to the selections of reasonable regression coefficients δ f wd i , i = 0, 1, 2. By comparing rational function-based regression with multiple selections of regression coefficients and linear regression, we also noticed that rational function-based regression produces higher correlation (considering PPMCC metric and contours) and lower regression errors (considering RMSE metric).
It should be noted that contours are intuitive visual indicators of correlation, while PPMCC and RMSE are statistical indicators of correlation or error. From Figures 11 and 13, Winkel Tripel projection has lower RMSE (about 1.61 • at best) between ω and ω f wd than RMSE for Collignon projection (about 3.21 • at best); however, ω f wd for Winkel Tripel projection has more visual deviations from the ω metric than deviations of the metrics for the Collignon projection.

Determination of Regression Coefficients
From the above comparison, there are different choices available for rational functionbased regression coefficients δ f wd i (see Table 3), i = 0, 1, 2. Furthermore, the proposed rational function-based regression is obviously insensitive to the selection of these reasonable regression coefficients. Now let us determine the rational function-based regression coefficients by considering conformal map projections whose ω = 0 or ρ d = 1.

Correlation and Error Analysis of Angular Distortion
Rational function-based regression (see Section 4.2) transforms the FWD-GCA-based metric ρ f wd (see Section 3) to ω f wd , an approximation to differential-based ω; therefore, we calculated ω, ω gc , and ω f wd to evaluate the effectiveness of the angular distortion metric ω f wd for all map projections we studied.
Tables S3-S6 in supplementary materials summarize PPMCC and RMSE between ω and ω gc (or ω f wd ) for more than 80 map projections that were investigated in this study.
In Table S3, cylindrical, conic, and azimuthal map projections used linear regression with coefficients β f wd 1 = 2 and β f wd 0 = −1. In Table S4, non-cylindrical, non-conic, and non-azimuthal map projections use rational function-based regression with coefficients The GCA-based metric cannot be calculated for map projections without inverse equations. As PROJ library does not provide inverse equations for all map projections studied, Table S5 provides the FWD-GCA metric only for those map projections whose inverse equations are not provided by PROJ library. Table S6 lists errors of angular distortion metrics for conformal map projections, where the maximum angular distortion is theoretically 0 • .
From Tables S3-S6, it can be inferred that both GCA-based and FWD-GCA metrics have fairly similar results to differential metric. Moreover, the proposed FWD-GCA metric appears to be more consistent with differential-based metrics than GCA-based metrics, although there is no significant difference for both metrics. The results also imply that the proposed FWD-GCA metric is an alternative to the GCA-based metric and that it can be used to calculate the appropriation of the maximum angular distortion metric. The worst statistical cases of the FWD-GCA metric in Tables S3-S6

Correlation and Error Analysis of Area Distortion
Area distortion is also significant for map projections. We calculated PPMCC and RMSE between area distortion metrics s and s gc (see Equation (22), or s f wd , see Equation (48)), and the results are listed in Tables S7-S9 in supplementary materials.
Metric s gc can be calculated for map projections with inverse equations; therefore, Table S7 provides the result for such types of map projections.
Metric s gc cannot be calculated for map projections without inverse equations; therefore, metric s f wd for these map projections was calculated and evaluated (see Table S8).
There are many equal area map projections; theoretically, s is equal to 1 for them. We calculated RMSE between metrics for these map projections (see Table S9).
Compared to the angular distortion, the proposed forward version of the area distortion s f wd (see Equation (48)) is highly correlated with the differential-based metric s (see Equation (9), or, precisely speaking, s f wd is almost the same as s).

Evaluation and Comparison of Map Projections
As the proposed metric and regression method proved to be valid for distortion evaluation for map projections, we used that metric to compare certain map projections. Figures 15 and 16 demonstrate this comparison. Due to the symmetry, only a quarter area (0 ≤ ϕ ≤ π/2 and 0 ≤ λ ≤ π) on the projection plane is plotted in Figures 15 and 16. The scalar bar of area distortion is plotted in the right side of Figure 15. Yellow represents no area distortion (s = 1), warm colors represent deflation of area (s < 1), and cool colors represent inflation of area (s > 1).
Certain map projections have the same or similar shapes, but their distortions are dissimilar. A comparison of the distortions of these similar map projections is beneficial in choosing the appropriate map projections for specific applications and controlling distortions under reasonable ranges.
As shown in Figure 15, the Hammer and Aitoff projections have the same shape (2:1 ellipses). From Figure 15a, we can observe that the Aitoff projection has a lower angular distortion than the Hammer projection in high latitude areas. However, the Hammer projection (equal area, not provided in Figure 15a) has a lower area distortion than the Aitoff projection.
From Figure 15b, we observe that the Apian Globular I and Ortelius oval map projections have similar distortions but the shapes of these two map projections are slightly different. In contrast, while Bacon Globular and Ortelius oval projections have very similar shapes, the distortions of these two map projections are obviously different.  It should also be noted that the PROJ library does not provide inverse equations for map projections that are listed in Tables S5 and S8. This means that it fails to calculate GCA-based metrics for these two map projections. In addition, calculating distributions and contours, as shown in Figures 15b and 16, is also inconvenient because one cannot directly establish approximately one million sample points on a projection plane, as shown in the map projections of Figures 9-14. As an alternative, we used the HEALPix algorithm to generate approximately three million sample points (12 × N 2 side , where N side = 512) on sphere S 2 , and calculated the metrics ρ d , ω, s, ρ f wd , ω f wd , and s f wd for these three million sample points. Then, we used cubic interpolation for the unstructured 2-D grid (i.e., on the two-dimensional projection plane) to calculate the metrics and contours (see Figures 15b and 16) for approximately one million sample points on the projection plane.

Discussion and Conclusions
We constructed and validated a forward version of the great circle arcs-based (or FWD-GCA-based) metrics ρ f wd and s f wd for evaluating angular and area distortion of map projections. The proposed metric is uniform irrespective of whether inverse equations for map projections exist or not, or the type of forward equations used (e.g., explicit or implicit). Although the differential-based metric also has unified formulae, the differential calculation may be complicated for certain map projections, and Newton-Raphson may be needed for some map projections with implicit equations. The use of ρ f wd avoids manual or auto calculation of the derivation, especially for map projections with singularity points of derivation, or implicit, interrupted, or non-differentiable forms of forward equations.
Compared to the differential-based metric and GCA-based metric, this study simplifies the process of calculating the distortions in map projections and expands the evaluation scope of map projections that can be evaluated by distortion metrics. The FWD-GCA-based metric is not dependent on inverse equations of map projections, and the only prerequisite for it is the existence of forward equations of the map projections. This implies that the proposed metric ρ f wd has a higher feasibility and usability than the differential metric ω.
Using rational function-based regression and linear regression, we could further estimate the maximum angular distortion ω f wd . Experimental results reveal that ω f wd , which is derived from ρ f wd , is a good approximation to ω. For non-cylindrical, nonconic, non-azimuthal, and non-conformal map projections studied in the experiments, the average errors of angular distortion were less than 1.778 • (i.e., the average of values in the last column in Tables S4 and S5 in  In addition, there are still certain noticeable errors between FWD-GCA-based metric ω f wd and differential-based metric ω in areas of lower distortions. However, for areas of higher distortions, which are of significant interest to researchers or mapmakers, ω f wd is more consistent with ω. A limitation of the proposed metric is the boundary problem on the projection plane. When λ approaches ±π, four points on the sphere (see Figure 3a) may be located on either side of the prime meridian for several map projections, and the corresponding four points on the projection plane (see p i in Figure 3b, where i = 1, 2, 3, 4) may be separated far away due to the projection. The above situation may incur errors in the process of vector interpolation (see Figures 3b,c). Therefore, ±π should be avoided for λ. In our implementation, approximately 200k sample points, generated by HEALPix, were shifted by 0.003 radian along the longitudinal direction. To calculate the metric of a sample point when λ = ±π for several map projections, the sample point can also be slightly shifted (e.g., by 1 × 10 −5 radian) along the longitudinal direction for an approximate calculation.
Additionally, the existing PROJ library provides coordinate transformation functions; and the proposed FWD-GCA metric employs PROJ to evaluate the distortion of map projections. However, the proposed FWD-GCA metric is not restricted to map projections provided in PROJ library; we can evaluate any map projections with forward equations by the FWD-GCA metric. As PROJ library provides convenience for the calculation of distortions, a practical implication of this study is the use of existing libraries. Apart from PROJ, there are several open or closed source softwares, such as NASA's G.Projector and Mapthematics LLC's Geocart. These mapping softwares support a greater number of map projections and provide transformation functions for raster images; however, they do not provide explicit mathematical functions for coordinate transformations. In the future, we will investigate practical numerical approaches to evaluate distortions based on raster images and further improve the differential-independent metric proposed in this study.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ijgi10100649/s1, Code S1: Direct calculation of differential-based metrics in Julia; Code S2: Iterative calculation of differential-based metrics in Julia Code S3: Julia code for FWD-GCAbased metrics; Table S1: Coefficients of linear (see Equation (50)) and rational function-based (see Equations (55)) regressions for each map projection; Table S2: Regression coefficients calculated by using a limited number of map projections; Table S3: Correlation (see Equation (58)) and errors (see Equation (59)) between differential-based metric ω and GCA-based metric ω gc (or FWD-GCAbased metric ω f wd ) for cylindrical, conic, and azimuthal map projections; Table S4: Correlation and errors between differential-based metric ω and GCA-based metric ω gc (or FWD-GCA-based metric ω f wd ) for map projections that have inverse equations; Table S5: Correlation and errors between differential-based metric ω and FWD-GCA-based metric ω f wd for map projections that do not have inverse equations; Table S6: Errors between differential-based metric ω and FWD-GCA-based metric ω f wd for conformal map projections; Table S7: Correlation and errors between differential-based and GCA-based area distortion metrics (or differential-based and FWD-GCA-based area distortion metrics) for map projections that have inverse equations; Table S8: Correlation and errors between differential-based and FWD-GCA-based area distortion metrics for map projections that do not have inverse equations; Table S9: Errors between differential-based and FWD-GCA-based area distortion metrics for equal area map projections.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Derivation and Implementation of Differential-Based Maximum Angular Distortion in Julia
Julia allows use of Unicode characters (such as ϕ, λ) in functions and variables; it also allows defining a mathematical function as a compact form of a callable function (e.g., x(ϕ, λ) = σ(ϕ) * sin( (ϕ, λ)) in Code S1 in supplementary materials). The above callable function can be used conveniently by ForwardDiff.jl, a Julia package, for calculating partial derivatives.
Due to the above reasons and the speed and flexibility offered by the language, we defined and used forward equations of map projections in Julia for further calculation of differential metrics, and also implemented the core algorithm of GCA-based and FWD-GCA-based (see Code S3) metrics in Julia.
Code S1 provides the Julia code for calculating differential metrics ω and s for the Bonne (ϕ 1 = 30 • ) projection. Most map projections could use the function 'calc_diff_generic' mentioned in Code S1 to calculate ω and s.
Code S2 provides the Julia code with Newton-Raphson iterations for calculating differential metrics ω and s for the Mollweide projection that do not have the explicit form of forward equations. Seventeen map projections with the implicit form of forward equations in our experiment, such as Boggs Eumorphic, Eckert IV, Eckert VI, Goode homolosine, Nell, and Putnin . š P 2 projections, used the function 'calc_diff_iterative' mentioned in Code S2 to calculate ω and s. Now, we give an explanation and derivation on map projections with implicit forward equations and Newton-Raphson iterations for solving implicit equations and obtaining the derivatives.
Equations (A1) and (A2) provide the forward equations of the Mollweide projection by using the prime meridian as the central meridian and using a unit sphere for simplicity. As the Mollweide projection does not have the explicit form of forward equations, we used the Newton-Raphson iterations to calculate the differential metric.