Plotting the map projection graticule involving discontinuities based on combined sampling

This article presents a new algorithm for interval plotting the projection graticule on the interval Ω = Ωφ × Ωλ based on the combined sampling technique. The proposed method synthesizes the uniform and adaptive sampling approaches and treats the discontinuities of the coordinate functions F,G. A full set of the projection constant values represented by the projection pole K = [φk, λk], two standard parallels φ1, φ2 and the central meridian shift λ0 are supported. In accordance with the discontinuity direction it utilizes a subdivision of the given latitude/longitude intervals Ωφ = [φ,φ], Ωλ = [λ, λ] to the set of disjoint subintervals Ω k,φ,Ω g k,λ forming tiles without the internal singularities, containing only “good” data; their parameters can be easily adjusted. Each graticule tile borders generated over Ω k = Ω g k,φ ×Ω g k,λ run along the singularities. For combined sampling with the given threshold α between the adjacent segments of the polygonal approximation, the recursive approach has been used; meridian/parallel offsets are ∆φ,∆λ. Finally, several tests of the proposed algorithms are involved.


Introduction
Maps are an essential part of our history and cultural heritage; a close attention is paid to their study and research.Working with the map content, its geometric and spatial characteristics described by the map projection cannot be ignored.The map projection P is defined by the coordinate functions F (ϕ, λ), G(ϕ, λ) of two independent variables, the latitude ϕ, and the longitude λ.F, G may have a different form; they may not be continuous, or their analytic form may not exist.To avoid the discontinuities, the interval Ω = Ω ϕ ×Ω λ will be divided into the disjoint set of k "good" subintervals Ω g k ; their "boundaries" run along the singularities.A current approach concentrated on uniform sampling of the meridians and parallels with the steps δϕ, δλ may not be sufficient.Despite its popularity, the equally spaced points cannot describe the meridian/parallel course without errors; the problems of undersampling or oversampling are common.Because of the complex/straight shapes of the meridians and parallels, this technique is not generally recommended.Adaptive sampling brings several benefits, it adapts to a different curvature of the function, reduces the amount of data and provides a natural and smooth plot of the function without jumps and unnatural breaks.This technique is popular in computer graphics; recall the deCasteljau or Chaikin's algorithms for the curve approximation.A combination of the uniform and adaptive sampling techniques synthesizes their advantages, which is discussed in [5].
Taking into account the facts mentioned above, this technique may provide a smooth and natural depiction of the graticule which is less data-redundant.A full set of the projection T. Bayer: Plotting the map projection graticule with discontinuities constant values represented by the projection pole K = [ϕ k , λ k ], two standard parallels ϕ 1 , ϕ 2 , and the central meridian shift λ 0 is supported.They have a strong influence on the shape of the graticule.The last parameter makes the process of the graticule reconstruction more difficult; an intersection of the meridian m(λ 0 ) with the parallel p(ϕ) as well as with the meridian m(λ) need to be determined.Another essential subproblem is represented by the automatic detection of the discontinuities; this issue refers to many coordinate functions F, G.
Finally, all algorithms will be tested on the real cartographic data represented by several projections.To illustrate the behavior, properties, and drawbacks of the proposed methods, the projections with more singularities will be preferred; the Fournier I. projection represents a typical candidate.Apart from current outcomes in the cartographic or GIS software tools, several specific applications may occur.Let us mention the problem of the unknown map projection analysis.Providing the visualization of the results of the detection algorithm illustrates its efficiency, which, from the "raw numbers", may not be apparent.This paper is organized as follows.In Section 3 the mathematical background of the solution is described.Section 4 is devoted to the proposed combined sampling algorithm, the properties of which in Section 5 will serve to analyze its behavior on real data.

Related Work
The polygonal approximation of the parametric curve based on adaptive sampling is mentioned in several papers: the refinement criteria in [12], the approximation by polygonal curves in [7], the spatial approximation of implicit curves in [17], [10], the affine arithmetic working in triangulated models in [21].The curve interpolation algorithms preserving its speed are described in [23], [25], [26], the adaptive-feed rate technique in [27].However, the map projections are never defined by the implicit equations.There are several papers focused on the approximation and drawing the functions or surfaces containing discontinuities.The interpolation and approximation of the piece-wise smooth functions are discussed in [2], [9].Several detection methods for the localization of the singularities of surfaces can be found in [13], [22], [19], [11], [18], [1], [8], the combined sampling technique in [5].The reconstructed projection graticule and its plot may be utilized in many cartographic or GIS applications.An interesting application is represented by the analysis of the unknown map projection described in [6], [4] providing the visualization of results of the detection algorithms, typically the reconstructed projection graticule.

Map projection singularities
During the graticule construction, several types of singularities occur; see Fig. 3.2.The sampling algorithm should be adapted to these facts.A meridian or a parallel may intersect a singularity or may be coincident.There are several simple strategies for handling and detecting the discontinuities; their overview can be found in [15], [3], [16], [5].Our approach is based on the LR criterion described in [20].
Coincidence with the infinite singularity.If the meridian m(λ) coincides with the infinite singularity c, λ = c, then, For Case B, the original meridian m(λ) is replaced with two new meridians m − (λ → c − ), and m + (λ → c + ).Analogously, Cases C, D occur, if a parallel p(ϕ) coincides with the infinite singularity c, Figure 3.1: Discontinuities of the coordinate functions F, G and asymptotes: the removable discontinuity (Hassler projection), the infinite discontinuity (gnomonic projection) at ϕ = 0, the infinite discontinuity (Apianus and Nicolosi projections) at λ = 0.
For Case C, the originally projected parallel p(c) may be replaced with the average for Case D the parallel p(ϕ) is replaced with two new parallels p − (ϕ → c − ), and p + (ϕ → c + ).
The removable discontinuities may be treated analogously to Cases A, C, the jump discontinuities analogously to Cases B, D. An example of the coincidence with the discontinuities can be found in Fig. 3.2.
Intersection of the singularity.The singularity may not affect the entire meridian or parallel, but only some inferior points.The meridian m(λ) intersects the infinite singularity c at the point representing Case F. For Case E, the projected meridian point q i is replaced with the average Treating Case F is more difficult, the entire meridian m(λ) needs be split into two parts at c T. Bayer: Plotting the map projection graticule with discontinuities projected separately.Analogously, Cases G, H occur, if a parallel p(ϕ) intersects the infinite singularity at the point lim For Case G, the projected parallel point q i is replaced with the average however, for Case H, the parallel p(ϕ) needs to be split into two parts at c The removable discontinuities may be treated analogously to Cases E, G, the jump discontinuities analogously to Cases F, H.

Detection of singularities.
For parametric functions F, G it is necessary to recognize whether a singularity refers to ϕ or λ coordinate.Therefore, each sampled point q i = [ϕ i , λ i ] and its boundary B(q i , ε) will be checked for a discontinuity.Suppose to be 5 sampled points.If any point of B ϕ i is singular, the singularity has a direction of the meridian m(λ i ) and analogously for B λ i .
For the oblique aspect of P, represent the projections of B(q i , ε) to the (ϕ , λ ) directions and q ϕ i+l = [ϕ i + lh, λ i ] ∈ B ϕ i , q λ i+l = [ϕ i , λ i + lh] ∈ B λ i are the 5 sampled points.If all points of B ϕ i are singular, the singularity has a direction of the meridian m(λ i ) and analogously for B λ i .While B ϕ i , B λ i refer to directions of the geographic meridian/parallel, B ϕ i , B λ i are aligned with the transformed meridian/parallel directions, see Fig. 3.3.
Depending on the direction, the coordinates x i , y i are labeled with ϕ, λ.The infinite discontinuities of coordinates x ϕ i−k , y ϕ i−k in the latitude/longitude directions are tested by comparing with the coordinate thresholds x, y.Subsequently, the jump discontinuity is detected using LR criteria then, F or G (or both) are probably not smooth at q i ; a discontinuity is aligned with the λ direction.Otherwise, if LR x (ϕ i + lh, λ i ) > LR ∨ LR y (ϕ i + lh, λ i ) > LR, a discontinuity is aligned with the ϕ direction, where LR = 0.8 represents the threshold.If a discontinuity is detected, the backward projection of B ϕ i , B λ i to B ϕ i , B λ i is performed; see Fig. 3.3.In practice, only the internal points q ϕ i+l of B ϕ i , and q λ i+l of B λ i are converted to the normal aspect using Eqs.3.4, 3.5 and their direction is checked.Subsequently, the meridian/parallel is shifted or split at T. Bayer: Plotting the map projection graticule with discontinuities O S

Intersection of the great circle and meridian/parallel arcs
Finding an intersection of the great circle arc and the meridian/parallel arc on the sphere represents an essential problem of the graticule construction algorithm.While p(ϕ), m(λ) refer to the normal aspect, commas in the labels indicate the oblique aspect.
The parallel p(ϕ c ) sampled on where λ 1 ≥ λ 2 , and ϕ 1 = ϕ 2 = ϕ c .The interval will be split into the three parts given by the longitude subintervals Depending on λ k , the longitude interval Ω λ will be split into three parts otherwise, Geoinformatics FCE CTU 17(2), 2018 T. Bayer: Plotting the map projection graticule with discontinuities For λ 0 = 0 the meridians m (0) and m(λ c ) have just two intersections at the poles: •] ≡ B; the latitude interval Ω ϕ of the meridian m(λ c ) does not need to be split, see Fig. 3.4.If λ 0 = 0, at most one intersection, M 3 = [ϕ 3 , λ c ] exists.To avoid the singularity, the meridian m(λ c ) sampled in Ω ϕ = [− π 2 , π 2 ] will be split into two parts given by the latitude subintervals Intersection of two planes.Suppose the planes ρ 1 (n 1 ), ρ 2 (n 2 ) given by the normal vectors which lies on both planes ρ 1 , ρ 2 will be chosen to minimize the norm [14] min It leads to the system of linear equations with the Lagrange multipliers λ 1 , λ 2 and the solution s = A −1 b, where Because A is sparse and symmetric, its inversion can be found in the analytic form.Let us put and the elements of The direction n of the line of intersection is given by the cross product Geoinformatics FCE CTU 17(2), 2018 T. Bayer: Plotting the map projection graticule with discontinuities Intersection of the line and sphere.It is necessary to test, whether the line of intersection intersects a sphere S 2 of the radius R with the center It leads to the quadratic equation for α, where For the discriminant D > 0, there are two intersections for D < 0 the line does not intersect the sphere.The Cartesian coordinates of intersections need to be converted to the spherical, where

The transformation to the geocentric coordinates
The spherical coordinates [ϕ i , λ i ] of points P 5 − P 13 are transformed to the Cartesian geocentric coordinates x i , y i , z i , where x i = cos ϕ i cos λ i , y i = cos ϕ i sin λ i , z i = sin ϕ i , i = 5, ..., 13.

Intersection of the line and sphere
Let S 2 be the sphere of the radius R = 1 centered at C = [0, 0, 0].Depending on the sign of the discriminant D, evaluate the intersections planes 2 and 3 .The planes 1 , 3 have two intersections denoted as

Transformation to spherical coordinates
Finally, the Cartesian coordinates of the intersections M 1 , M 2 , M 3 will be converted to the spherical.

Graticule reconstruction using the combined sampling
Let and Ω λ = [λ, λ], be the subdomain, inside which the polygonal approximation of the meridians/parallels is constructed.The points ), and analogously for the parallel.The combined sampling procedure represents a generalization of the curve sampling procedure discussed in [5].It starts with the uniform sampling procedure, which after d steps transforms to adaptive sampling.A refinement criterion is based on the angular difference α i of segments formed by three adjacent points (p i−1 , p i , p i+1 ).

Graticule construction with singularities
The proposed algorithm based on the stack S implementation is described in Alg. 1.The basic idea is to find a set of disjoint subsets Ω g ϕ , Ω g ϕ ⊆ Ω ϕ , and Ω g λ , Ω g λ ⊆ Ω λ , where Ω g = Ω g ϕ × Ω g λ , containing "good" data without singularities that allow for adaptive sampling.Over each Ω g j , the borders of which run along discontinuities the graticule fragment is constructed.The entire graticule is put together from these fragments.
The point q i = [ϕ i , λ i ] is "good" if no singularity at F (q i ) and G(q i ) occurs.Unlike the one-dimensional version of the problem, Ω is partitioned in two orthogonal directions (ϕ, λ).Let the j − th interval Ω j = Ω ϕ,j × Ω λ,j , where Ω ϕ,j = [ϕ j , ϕ j ], Ω λ,j = [λ j , λ j ], contain a singularity c, c ∈ Ω, with the unknown direction stored in the stack S, and ε, ε > 0, be the numerical threshold.Suppose ∆ϕ to be the latitude offset between parallels and ∆λ to be the longitude offset between meridians.

T. Bayer: Plotting the map projection graticule with discontinuities
Algorithm 1 Combined sampling with the singularities, the stack implementation.
while S = ∅ do 5: if (s > s) then 21: All adaptively sampled meridian/parallel points q i are checked for the discontinuities.If a discontinuity c is found, it is classified, and its direction is determined.Depending on c value, the lower/upper bound of Ω ϕ,j or Ω λ,j is shifted, or, a split to two disjoint intervals is performed: Otherwise, the polygonal approximation L λ,j of the meridian m(λ) or L ϕ,j of the parallel p(ϕ) is constructed.All disjoint polygonal approximations L j are stored in the lists of meridians L λ and parallels L ϕ .Recall s to be the amount of Ω ϕ,j , and Ω λ,j splits.The procedure can be summarized as follows: 1.

The initial phase
Initialize the empty stack S = ∅, Ω g = Ω and push S ← Ω. Set the number of splits to s = 0.

T. Bayer: Plotting the map projection graticule with discontinuities
Algorithm 2 The construction of meridians.

Recursive steps
Repeat the following steps until S is empty: (a) Pop the actual interval Ω j Pop the actual good interval Ω j ← S from S, where Ω g j = Ω ϕ,j × Ω λ,j , and

(b) Create the empty lists
Create new empty lists L λ,j = ∅, L ϕ,j = ∅ storing the polygonal approximation of the meridians and parallels.

(c) Create the intervals for feasible splitting
Initialize the newly created intervals Ω j,1 = Ω ϕ,j,1 ×Ω λ,j,1 , and and ] used for feasible splitting of Ω g j .Set the amount of newly created intervals as k = 0.
(d) Create the temporary polygonal approximation L ϕ,j , L λ,j of the graticule Create the temporary polygonal approximation of the meridians and parallels on Ω g j using adaptive sampling stored in L ϕ,j , L λ,j .

(e) Copy the temporary polygonal approximation
If no discontinuity appears, add L λ,j to L λ : L λ ← L λ,j , and L ϕ,j to L ϕ : L ϕ ← L ϕ,j , and go to step (a).Otherwise, c represents the discontinuity of the given type and direction detected.The following cases, when c = ϕ i , or, c = λ i , must be treated in Steps i-iii).
iii.If c = ϕ i ∨ c = λ i , do the following steps.If s > s, the maximum allowed recursion depth is exceeded without a reasonable solution, clear the polygonal approximations L ϕ , L λ .Otherwise, if k > 0, at least one new interval needs to be created; push Ω j,1 to the stack: S ← Ω j,1 .If k > 1, push the second interval Ω j,2 to the stack: S ← Ω j,2 .
For the implementation, see Alg. 1.Its design is robust against common numerical failures.The singularity value c is stored in the thrown exceptions SingularityLatError, SingularityLonError, which are derived from the base class SingularityError.The splitting procedure is realized in accordance with the singularity direction.

Generate meridians.
The procedure creates all meridians inside the given interval Initially, the bounding meridians m(λ j ), and m(λ j ) are created.To find the longitude λ start of the first meridian m(λ start ), which is the smallest multiplier of ∆λ higher than λ j , and the last meridian m(λ end ), which is the largest multiplier of ∆λ lower than λ j , use Other meridians m(λ j,i ) are generated inside the interval [λ start , λ end ] so that λ j = λ j,start + i∆λ, λ start < λ < λ end .Due to possible discontinuities, the meridian consists of several parts.
If a discontinuity at c = ϕ j is found, the meridian needs to be split at c.If a discontinuity at c = λ j is detected, a meridian needs to be shifted.This is achieved by splitting Ω j at c = ϕ j , or, at c = λ j with the additional shift ε.Otherwise, a current meridian m(λ j ) is constructed.Finally, the meridians passing along the pole m(λ k ± ε) and the opposite meridians m(λ k ± π ∓ ε), forλ k ≶ 0 are added.The procedure is summarized in Alg. 2.

Generate parallels.
The procedure generates all parallels inside the interval Initially, the bounding parallels p(ϕ j ), and p(ϕ j ) are created.The latitude ϕ start of the first parallel p(ϕ start ) is the smallest multiplier of ∆ϕ higher than ϕ j , the last parallel p(ϕ end ) is the largest multiplier of ∆ϕ lower than ϕ j .Other parallels p(ϕ j,i ) are generated inside the interval [ϕ start , ϕ end ] so that ϕ j = ϕ start + i∆ϕ, ϕ start < ϕ < ϕ end .If a discontinuity at c = λ j is found, the parallel needs to be split at c.If a discontinuity at c = ϕ j is detected, the parallel needs to be shifted.This is achieved by splitting Ω j at c = ϕ j , or, at c = λ j .Otherwise, a current parallel p(ϕ j ) is constructed.The procedure is summarized in Alg. 3.

Generate a meridian.
Combined sampling provides a polynomial approximation of the meridian m(λ) by the refinement α and the recursion depth d over the interval Ω ϕ,j = [ϕ j , ϕ j ].
T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 3 The construction of parallels.
If an intersection M 3 is found, the splitting procedure is undertaken; see Fig. 3.4 and Alg. 4.

A combined sampling of the fragments.
The sampling methods will be discussed in Sec.4.2.
Generate a parallel.The parallel p(ϕ) is intersected by the meridian m (λ 0 ) in at most two points.The original interval Ω λ,j needs to be split along the intersection points ; the solution is found from the intersection of two planes (see Sec. 3.2).In general, the following cases are solved: For practical computation, the coordinates are sorted so that λ 1 ≤ λ 2 ; see Alg. 5.

Algorithm 4
The construction of a meridian from the fragments.

Combined sampling of meridian and parallel fragments
Three sampling techniques denoted as S1, S2, S3 representing a natural extension of the sampling methods described in [5] will be illustrated on the meridian/parallel approximation.

Sampling method S1.
It provides a polynomial approximation of the meridian m(λ) by the refinement criteria α and the recursion depth d.Both Ω ϕ,j,1 , Ω ϕ,j,2 are partitioned into four disjoint subintervals with the randomly shifting borders.

The initial phase
Let L λ,j,k = ∅ be the empty set.Compute Set the recursion depth d = 1.The procedure is summarized in Alg. 6.

The recursive procedure
Enter the recursive procedure and do the following substeps: into the approximate quarters where r 1 , r 2 , r 3 are the random numbers inside the interval [0.45, 0.55].

Algorithm 5
The construction of a parallel from the fragments.

Add the last point
Add the last point p b to the polynomial approximation of m(λ) : L λ,j,k ← p b and finish the adaptive sampling procedure.
The procedure is summarized in Alg. 7. Unfortunately, this technique does not involve the angles α between the adjacent segments from the previous iteration.While the newly created segments of the polygonal approximation of m(λ) can fulfill the condition of α i ≤ α, their joining to the previous/next segments may not hold this requirement.In general, this issue affects curves of the complex shape, where the average maximum value α i .= 1.25α.

Sampling method S2.
The first improvement is represented by the refinement of the recursive condition.The original conditions and α 3 > α, are replaced with the new identical for all recursive calls.This approach brings a slightly increased amount of the recursive calls, while the average maximum value falls to α i .= 1.15α.In most situations, this issue does not significantly affect the smoothness of the polynomial approximation and its perception by the user.However, exceeding the threshold may not always be acceptable.

Sampling method S3
. Unlike S1, S2, the improved sampling method S3 takes into account angles between the newly created and adjacent segments from the previous iteration.

T. Bayer: Plotting the map projection graticule with discontinuities
Algorithm 7 Combined sampling of the meridian points, the recursive phase, method S1.P, L λ,j,k , a, b, x a , y a , x b , y b , d, d, d, , α)) throw LatSingularityException (ϕ i ), i = 1, 2, 3 throw LonSingularityException (λ), i = 1, 2, 3 10: : : L λ,j,k ← P oint(x 1 , y 1 ) L λ,j,k ← P oint(x 3 , y 3 ) The approximate quarters of Ω ϕ,j,k denoted as Ω ϕ,j,k,l , l = 1, ..., 4, have the form of The idea of modification is straightforward, an information about previous and subsequent quarters of Ω ϕ,j,k,l is stored While the first quarter lacks its predecessor, the last quarter lacks the successor.For a current recursive depth d, the segment (p a,k , p b,k ) of the polygonal approximation is split to the quarters, where p in the recursive depth d + 1, but also between the adjacent segments where p 3,k−1 is the projected lower bound of Ω ϕ,j,k−1,4 (a projected third quarter of the interval Ω ϕ,j,k−1 previous to Ω ϕ,j,k ) with the latitude p 1,k+1 is the projected first-quarter Ω ϕ,j,k+1,1 of the next interval Ω ϕ,j,k+1 with the latitude and r 0 , r 4 are the random numbers generated according to the principles mentioned above.The modified condition for the recursive call has the form of In general, the improved performance of S3 occurs only in specific situations when the sampled function has a complex shape (meridians of Bonne projections).However, the significant disadvantage is represented by the increased amount of sampled points in the highly curved areas.The recursive step has a slightly different form: stop the recursive procedure and go to Step 3.

For a given
, and successor Ω ϕ,j,k+1 = [a k+1 , b k+1 ], the interval is split by three points, the positions are close to the quarters of the interval Geoinformatics FCE CTU 17(2), 2018 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 8 Combined sampling of the meridian points, the recursive phase, method S3.
, λ] occurs, throw a new exception according to the discontinuity direction.
(b) Otherwise, evaluate the function values T. Bayer: Plotting the map projection graticule with discontinuities The discontinuity treatment without its classification using the proposed algorithm in polyconic projection; the doubled equator is visible under magnification.
), and the recursive depth d.When a meridian is not sufficiently smooth, or d ≤ d, it needs to be refined; the recursive subdivision is necessary.
(d) If the condition given by Eq. 4.2 is held, call the recursive procedure with the increased depth (e) Add the new point p 1,k to the polynomial approximation of m(λ) : L λ,j,k ← p 1,k .
(f) If the condition given by Eq. 4.2 is held, call the recursive procedure with the increased depth (g) Add the new point p 2,k to the polynomial approximation of m(λ) : L λ,j,k ← p 2,k .
(h) If the condition given by Eq. 4.2 is held, call the recursive procedure with the increased depth (i) Add the new point p 3,k to the polynomial approximation of m(λ) : L λ,j,k ← p 3,k .
(j) If the condition given by Eq. 4.2 is held, call the recursive procedure for the last interval The procedure is summarized in Alg. 8.
A combined sampling of the parallel.Analogously, the parallel p(ϕ) is sampled over the interval Ω λ,j = [λ j , λ j ].The polygonal approximation consists of two steps: a partition of p(ϕ) to the fragments and combined sampling of the fragments.
Each parallel fragment needs to be sampled, Ω ϕ,j,1 , Ω ϕ,j,2 , Ω ϕ,j,3 are partitioned into four T. Bayer: Plotting the map projection graticule with discontinuities disjoint subintervals.While for adaptive meridian m(λ) sampling are: λ = const, ϕ ∈ Ω ϕ,j , for parallel p(ϕ) sampling are: ϕ = const, λ ∈ Ω λ,j .Therefore, the randomly generated parallel points split the interval Ω λ,j,k to the approximate quarters The remaining steps are analogous.This method can treat all kinds of discontinuities, without specific knowledge of ε; for practical computations is ε = 0.001 • .It has a disadvantage -the duplication of lines along singularities if the discontinuities are not classified.For the discontinuity c = ϕ = 0 • in the polyconic projection, this issue is illustrated in Fig. 4.3.A hemisphere constructed in the normal aspect is composed of two sub-intervals [−π/2, 0), (0, π/2].Hence, the equator is doubled, which does not look aesthetically pleasing, but it is visible only under magnification.

Experiments and results
The principles mentioned above will be illustrated on the several tests.While the first analysis compares the behavior and properties of combined sampling, the second test reconstructs the Fournier I. projection graticule in the oblique aspect (several discontinuities in the coordinate functions are involved).Finally, the last experiment illustrates the splitting procedure efficiency depending on the threshold ε.The results are summarized in Tabs.1-5.

Comparison of uniform and combined sampling
During this test, the uniform and combined sampling techniques will be compared regarding the data representation compactness measured by the amount of the sampled meridian points n mer and parallel points n par .Additionally, the maximum angles α m , α p between the sampled meridian and parallel segments together with their mean values α m , α p are measured; these indicators depend on the sampling density.For the uniform sampling, the sampling steps δ ϕ , δ λ are fixed, while the adaptive sampling is driven by the sampling angle α.Our experiments will be carried out for combined sampling with α = 1 • , 2 • , 5 • , the recursive depth is d = 1, and for the uniform sampling where δ ϕ , δ λ will be set providing the analogous amount of the sampled points.8 projections are involved in testing: equidistant cylindrical, conic, azimuthal, Bonne and Hassler (ϕ 1 = 50 • ), Nicolosi, Littrow, and Adams I World conformal.The results are summarized in Tab. 1.
While S1, S2 exceed the α m , α p values (+20% for S1, +15% for S2), S3 provides values satisfying the given criteria (-7%) without a significant increment of sampled points.Cylindrical, conic and azimuthal projections illustrate the adaptive sampling efficiency; for the straight meridians/parallels, only 2% of points are sufficient for the shape estimation (uniform sampling cannot achieve this efficiency).
Concerning the criteria mentioned above the results are compared in Tab. 2. The sampling steps δ ϕ , δ λ have been set so that the amount of sampled points n mer , n par is analogous.In T. Bayer: Plotting the map projection graticule with discontinuities  the higher-curvature regions, the combined sampling technique provides a smoother approximation of the meridians/parallels.A typical example is represented by the Bonne, Littrow or Adams I. projections, where uniform sampling does not provide good results.The maximum angles between the sampled elements are more than five times higher.If meridians or parallels are represented by the circular arcs, the combined sampling does not bring any advantage and may cause a minor deterioration of the average values α m , α p (Nicolosi projection).Raise the steps δϕ, δλ reduces the advantage of adaptive sampling (and thus α) even for the highcurvature regions (see Bonne projection in Fig. 5.1).In general, uniform sampling preserves the curvature worse and provides more redundant data.

Detection of discontinuities: Fournier I. projection
The Fournier I. projection equations have a complex form and contain several discontinuities both in the latitudinal and longitudinal directions.In accordance with [24], the projection equations are written as follows.

Measuring the quantitative parameters of the algorithm
The algorithms have been tested on the construction of the several graticules with the shifted central meridian λ 0 .In all cases, the entire planisphere is covered by the graticule with the offset of ∆ϕ = ∆λ = 10 • .Only the map projections with the singularities are involved in testing: Fournier I., Van der Grinten I., Ortelius Oval, Nicolosi, Hassler.While the Hassler projection has the only singularity at ϕ = 0 • , the Ortelius projection is singular at λ − λ 0 = 0, and the Fournier I., Van der Grinten I., and Nicolosi projections add singularities at λ − λ 0 = ± π 2 .The graphical issues, amount of splits n s , processed tiles n t , interval resizing n r , deleted intervals n d , and possible failures of the algorithm will be investigated.

M 1 Figure 3 . 4 :
Figure 3.4: The intersections of the meridian m (λ 0 ), parallel p(ϕ c ), and the meridian m(λ c ) lead to the ϕ/λ interval partition into two/three subintervals; A represents the North Pole, K the transformed pole.

( a )
If d > d or b − a < ε, stop the recursive procedure and go to Step 3. (b) For a given Ω ϕ,j,k = [a, b], the interval is split by three points
call the recursive procedure with the increased depth d = d + 1 for the interval [ϕ 2 , x 3 ].(k) Add the new point p 3 to the polynomial approximation of m(λ) : L λ,j,k ← p 3 .(l) If α 3 > α, call the recursive procedure for the interval [ϕ 3 , b].

Figure 4 . 2 :
Figure 4.2: Illustration of the combined sampling method S3 for the meridian m(λ). α k and the supplementary vertices p 3,k−1 , p 1,k+1 .(c) For d ≤ d, this step begins with uniform sampling of the meridian points.If d > d, it transforms to the adaptive method.Check the refinement criteria

Figure 4 . 3 :
Figure 4.3: The discontinuity treatment without its classification using the proposed algorithm in polyconic projection; the doubled equator is visible under magnification.

Figure 5 . 2 :
Figure 5.2: Fournier I. projection with the highlighted singularities leading to the subdivision into the tiles; normal aspect of the projection.

Figure 5 . 3 :Figure 5 . 4 :
Figure 5.3: Reconstruction of the graticule from the tiles generated along the set of Ω g , Fournier I projection.

Table 2 :
Uniform sampling of the graticule with the determined steps δϕ, δλ compared to S3; the amounts of sampled points n mer , n par are preserved.
T. Bayer: Plotting the map projection graticule with discontinuities