Plotting the map projection graticule involving discontinuities based on combined sampling Plotting the map projection graticule involving discontinuities based on combined sampling Tomáš Bayer Department of Applied Geoinformatics and Cartography, Faculty of Science, Charles University, Czech Republic bayertom@natur.cuni.cz Abstract. 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 discontinu- ities 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 subdivi- sion of the given latitude/longitude intervals Ωϕ = [ϕ,ϕ], Ωλ = [λ, λ] to the set of disjoint subintervals Ωg 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 Ωg k = Ωg k,ϕ ×Ω g k,λ run along the singularities. For combined sampling with the given thresh- old α 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. Keywords: digital cartography; mathematical cartography; adaptive sampling; graticule; meridians; parallels; recursive approach; map projection; great circle; discontinuity; visual- ization; sphere. 1. 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 Geoinformatics FCE CTU 17(2), 2018, doi:10.14311/gi.17.2.3 31 http://orcid.org/0000-0001-6279-1926 http://dx.doi.org/10.14311/gi.17.2.3 http://creativecommons.org/licenses/by/4.0/ 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. 2. Related Work The polygonal approximation of the parametric curve based on adaptive sampling is men- tioned 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 de- scribed in [6], [4] providing the visualization of results of the detection algorithms, typically the reconstructed projection graticule. 3. Map projection and its properties Let S2 be a sphere of the radius R in R3 (reference surface, Earth), A = (0, 0, R), B = (0, 0,−R) North and South Poles of S2, and σ a plane. For a current point Q = [ϕ, λ] ∈ S2, different from A,B, and its image P ′ = [X,Y ] ∈ σ, the map projection P : S2 − {A,B} → σ, P(Q) = P ′, is defined by the coordinate functions F,G, of two independent variables ϕ, λ X = F (ϕ, λ), Y = G(ϕ, λ), (3.1) which are continuous with their first order partial derivatives, finite, and independent. A transformation between the normal and oblique aspects is performed using the laws of spher- Geoinformatics FCE CTU 17(2), 2018 32 T. Bayer: Plotting the map projection graticule with discontinuities ical trigonometry sinϕ′ = sinϕk sinϕ+ cosϕk cosϕ cos ∆λ, (3.2) tanλ′ = cosϕ sin ∆λ cosϕ sinϕk cos ∆λ− sinϕ cosϕk , (3.3) where ϕ′, λ′ are the spherical coordinates related to K. The inverse transformation has the form of sinϕ = sinϕk sinϕ′ − cosϕk cosϕ′ cosλ′, (3.4) tan ∆λ = cosϕ′ sinλ′ cosϕ′ sinϕk cosλ′ + sinϕ′ cosϕk , (3.5) where ∆λ = λ−λk. A meridian of the sphere S2 is a curve m(ϕ, λc), ϕ ∈ [−π 2 , π 2 ], λc = const; its image in P is m̃(ϕ, λc) = (F (ϕ, λc), G(ϕ, λc)). A parallel of the sphere S2 in R3 is a curve p(ϕc, λ), λ ∈ [−π, π], ϕc = const; its image in P is p̃(ϕc, λ) = (F (ϕc, λ), G(ϕc, λ)). In the simplified notation m(λ) represents a meridian of the longitude λ, p(ϕ) a parallel of the latitude ϕ and m̃(λ), p̃(ϕ) their projected variants. 3.1. 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, lim λ→c± F (ϕ, λ) =∞∨ lim λ→c± F (ϕ, λ) = −∞∨ lim λ→c± G(ϕ, λ) =∞∨ lim λ→c± G(ϕ, λ) = −∞, (Case A), or, lim λ→c± F (ϕ, λ) = ±∞∨ lim λ→c± F (ϕ, λ) = ∓∞∨ lim λ→c± G(ϕ, λ) = ±∞∨ lim λ→c± G(ϕ, λ) = ∓∞, (Case B) for any ϕ ∈ (−π/2, π/2). Using the left- and right-boundaries (c+ ε) and (c− ε) in Case A, the original projected meridian m̃(c) is replaced with m̃(c) = 0.5 [m̃(c− ε) + m̃(c+ ε)] . 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, lim ϕ→c± F (ϕ, λ) =∞∨ lim ϕ→c± F (ϕ, λ) = −∞∨ lim ϕ→c± G(ϕ, λ) =∞∨ lim ϕ→c± G(ϕ, λ) = −∞, or, if, lim ϕ→c± F (ϕ, λ) = ±∞∨ lim ϕ→c± F (ϕ, λ) = ∓∞∨ lim ϕ→c± G(ϕ, λ) = ±∞∨ lim ϕ→c± G(ϕ, λ) = ∓∞. Geoinformatics FCE CTU 17(2), 2018 33 T. Bayer: Plotting the map projection graticule with discontinuities −90 −60 −30 0 30 60 90 −10 −5 0 5 10 Hassler projection −90 −60 −30 0 30 60 90 0 −50 100 150 Gnomonic projection −180 −120 −60 0 60 120 180 0 50 100 150 Apianus projection F (λ ) [k m ] −180 −120 −60 0 60 120 180 −150 −100 −50 0 50 100 150 Nicolosi projection λ [deg]λ [deg] ϕ [deg] ϕ [deg] F (λ ) [k m ] F (ϕ ) [k m ] G (ϕ ) [k m ] a) b) c) d) 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 p̃(c) = 0.5 [p̃(c− ε) + p̃(c+ ε)] , (3.6) 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 discon- tinuities 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 qi = [ϕi = c, λ], if lim ϕi→c± F (ϕi, λ) =∞∨ lim ϕi→c+ F (ϕi, λ) = −∞∨ lim ϕi→c± G(ϕi, λ) =∞∨ lim ϕi→c± G(ϕi, λ) = −∞, (Case E), or, if lim ϕi→c± F (ϕi, λ) = ±∞∨ lim ϕi→c± F (ϕi, λ) = ∓∞∨ lim ϕi→c± G(ϕi, λ) = ±∞∨ lim ϕi→c± G(ϕi, λ) = ∓∞, representing Case F. For Case E, the projected meridian point qi is replaced with the average Xi = 0.5 (F (c− ε, λ) + F (c+ ε, λ)) , Yi = 0.5 (G(c− ε, λ) +G(c+ ε, λ)) . Treating Case F is more difficult, the entire meridian m(λ) needs be split into two parts at c m(λ) = 〈〈 m−(ϕi, λ), ϕ ≤ ϕi < c 〉 , 〈 m+(ϕi, λ), c < ϕi ≤ ϕ 〉〉 , Geoinformatics FCE CTU 17(2), 2018 34 T. Bayer: Plotting the map projection graticule with discontinuities A K m(l k )ºm (0) m (l 0 ) l k l 0 m(0) l k ±p m(l k ±p) m(±p) Figure 3.2: Coincidence with the singularities: the meridian m(λk) passing the transformed pole K, the meridian m(λk ± π) “opposite” the transformed pole, and the meridian m(±π); the azimuthal projection. projected separately. Analogously, Cases G, H occur, if a parallel p(ϕ) intersects the infinite singularity at the point qi = [ϕ, c = λi], lim λi→c± F (ϕ, λi) =∞∨ lim λi→c± F (ϕ, λi) = −∞∨ lim λi→c± G(ϕ, λi) =∞∨ lim λi→c± G(ϕ, λi) = −∞, or, if, lim λi→c± F (ϕ, λi) = ±∞∨ lim λi→c± F (ϕ, λi) = ∓∞∨ lim λi→c± G(ϕ, λi) = ±∞∨ lim λi→c± G(ϕ, λi) = ∓∞. For Case G, the projected parallel point qi is replaced with the average Xi = 0.5 (F (ϕ, c− ε) + F (ϕ, c+ ε)) , Yi = 0.5 (G(ϕ, c− ε) +G(ϕ, c+ ε)) , however, for Case H, the parallel p(ϕ) needs to be split into two parts at c p(ϕ) = 〈〈 p−(ϕ, λi), λ ≤ λi < c 〉 , 〈 p+(ϕ, λi), c < λi ≤ λ 〉〉 . The removable discontinuities may be treated analogously to Cases E, G, the jump disconti- nuities analogously to Cases F, H. Geoinformatics FCE CTU 17(2), 2018 35 T. Bayer: Plotting the map projection graticule with discontinuities q i q i-1 q i-2 q i+1 q i+2 q i+2 q i+1 q i q i-1 q i-2 p(φ) m(λ) m(λ‘) p(φ‘) m(λ) m(λ‘) p(φ) p(φ‘) B φ‘ B l‘i i Figure 3.3: Boundaries Bϕ′i , Bλ′i and their backward projections to the (ϕ, λ) directions. Detection of singularities. For parametric functions F,G it is necessary to recognize whether a singularity refers to ϕ or λ coordinate. Therefore, each sampled point qi = [ϕi, λi] and its boundary B(qi, ε) will be checked for a discontinuity. Suppose Bϕi = [ϕi − ε, ϕi + ε], Bλi = [λi − ε, λi + ε] to be the projections of B(qi, ε) to the (ϕ, λ) directions and qϕi+l = [ϕi + lh, λi] ∈ Bϕi , qλi+l = [ϕi, λi + lh] ∈ Bλi , where h = ε/2, l = −2,−1, 0, 1, 2, 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, Bϕ′i = [ϕ′i−ε, ϕ′i+ε], Bλ′i = [λ′i−ε, λ′i+ε] represent the projections of B(qi, ε) 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 xi, yi are labeled with ϕ, λ. The infinite disconti- nuities 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 LRx(ϕ′i, λ′i) = ∣∣F 2 r − F 2 l ∣∣ F 2 r + F 2 l , LRy(ϕ′i, λ′i) = ∣∣G2 r −G2 l ∣∣ G2 r +G2 l , where Fr = 3Fi− 4Fi+1 +Fi+2, Fl = 3Fi− 4Fi−1 +Fi−2, and Gr = 3Gi− 4Gi+1 +Gi+2, Gl = 3Gi − 4Gi−1 + Gi−2. If LRx(ϕ′i, λ′i + lh) > LR ∨ LRy(ϕ′i, λ′i + lh) > LR, l = −2,−1, 0, 1, 2, then, F or G (or both) are probably not smooth at qi; a discontinuity is aligned with the λ′ direction. Otherwise, if LRx(ϕ′i + lh, λ′i) > LR ∨ LRy(ϕ′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 qi = [ϕi, λi]. Geoinformatics FCE CTU 17(2), 2018 36 T. Bayer: Plotting the map projection graticule with discontinuities A KM 2 M 3 m(l c ) p(j c ) p(j 3 ) m (0) m (l 0 ) l k l c l 0 m(0) l 1 l 2 m(l 2 ) m(l 1 ) l=p M 1 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. 3.2. 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 Ωλ = [−π, π] intersects m′(λ′0) in at most two points M1 = [ϕ1, λ1], M2 = [ϕ2, λ2] where λ1 ≥ λ2, and ϕ1 = ϕ2 = ϕc. The interval will be split into the three parts given by the longitude subintervals Ωλ,1 = [−π, λ2), Ωλ,2 = (λ2, λ1), Ωλ,3 = (λ1, π]. If λ′0 = 0, then M1 = [ϕc, λk], M2 = [ϕc, λk ± π]. Depending on λk, the longitude interval Ωλ will be split into three parts Ωλ,1, Ωλ,2, Ωλ,3. If λk ≥ 0 Ωλ,1 = [−π, λk − π), Ωλ,2 = (λk − π, λk), Ωλ,3 = (λk, π], otherwise, Ωλ,1 = [−π, λk), Ωλ,2 = (λk, λk + π), Ωλ,3 = (λk + π, π]. Geoinformatics FCE CTU 17(2), 2018 37 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: M3 = [π2 , ·] ≡ A, M4 = [−π 2 , ·] ≡ B; the latitude interval Ωϕ of the meridian m(λc) does not need to be split, see Fig. 3.4. If λ′0 6= 0, at most one intersection, M3 = [ϕ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 Ωϕ,1 = [−π2 , ϕ3), Ωϕ,1 = (ϕ3, π 2 ]. Intersection of two planes. Suppose the planes ρ1(n1), ρ2(n2) given by the normal vectors n1 = (n1,x, n1,yn1,z), n2 = (n2,x, n2,yn2,z), the points Q1 ∈ ρ1, Q2 ∈ ρ2 and an arbitrary point Q0 = [x0, y0, y0]. The intersection point M = [xM , yM , zM ] which lies on both planes ρ1, ρ2 will be chosen to minimize the norm [14] min ‖M −Q0‖22 , subject to (M −Q1)n1 = 0, (M −Q2)n2 = 0. It leads to the system of linear equations with the Lagrange multipliers λ1, λ2 and the solution s = A−1b, where s =  xM yM zM λ1 λ2  , A =  2 0 0 n1,x n2,x 0 2 0 n1,y n2,y 0 0 2 n1,z n2,z n1,x n1,y n1,z 0 0 n2,x n2,y n2,z 0 0  , b =  2x0 2y0 2z0 x1n1,x + y1n1,y + z1n1,z x2n2,x + y2n2,y + z2n2,z  . Because A is sparse and symmetric, its inversion can be found in the analytic form. Let us put k0 = n2 1,x, k1 = n2 1,y, k2 = n2 1,z, k3 = n2 2,x, k4 = n2 2,y, k5 = n2 2,z, k6 = k0 + k1 k7 = k3 + k4, k8 = k3 + k5, k9 = k4 + k5, k10 = n1,xn2,x, k11 = n1,xn2,x, k12 = n1,zn2,z, k13 = n1,xn2,y, k14 = n1,xn2,z, k15 = n1,yn2,x, k16 = n1,yn2,z, k17 = n1,zn2,y, k18 = n1,zn2,x, k19 = k15n2,y, k20 = k17n2,z, k21 = k15n1,x, k22 = k13n1,y, k23 = k0 + k1 + k2, k24 = k3 + k4 + k5, then det(A) = 2(k0k9 − 2k10k12 + k2k7 − 2k11(k10 + k12) + k1k8), and the elements of A−1 are a11 = (k18 − k16)(k18 − k16) det(A) , a12 = (k14 − k17)(k18 − k16) det(A) , a13 = (k13 − k15)(k16 − k18) det(A) , a14 = 2(−k19 − k20 + n1,xk9) det(A) , a15 = 2[n1,yk15 − k22 + n1,z(k17 − k14)] det(A) , a22 = (k17 − k14)(k17 − k14) det(A) , a23 = (k13 − k15)(k17 − k14) det(A) , a24 = 2[n1,yk8 − n2,y(k10 + k12)] det(A) , a25 = 2[n1,xk13 − k21 + n1,z(k18 − k16)] det(A) , a33 = (k15 − k13)(k15 − k13) det(A) , a34 = 2[n1,zk7 − n2,z(k10 + k11)] det(A) , a35 = 2[n2,zk6 − n1,z(k10 + k11)] det(A) , a44 = −4k24 det(A) , a45 = 4(k10 + k11 + k12) det(A) , a55 = −4k23 det(A) . The direction n of the line of intersection is given by the cross product n = n1 × n2, its parametric equation is (x, y, z) = (xM , yM , zM ) + α(nx, ny, nz), where α ∈ [0, 1]. Geoinformatics FCE CTU 17(2), 2018 38 T. Bayer: Plotting the map projection graticule with discontinuities Intersection of the line and sphere. It is necessary to test, whether the line of in- tersection intersects a sphere S2 of the radius R with the center C = [xc, yc, zc], given by (xM + αnx − xc)2 + (yM + αny − yc)2 + (zM + αnz − zc)2 = R2. It leads to the quadratic equation for α, where a = n2 x + n2 y + n2 z, b = 2(xMnx − nxxc + yMny − nyyc + zMnz − zMzc), c = (xm − xc)2 + (ym − yc)2 + (zm − zc)2 −R2. For the discriminantD > 0, there are two intersectionsM1 = [x1, y1, z1], andM2 = [x2, y2, z2], determined from (x1,2, y1,2, z1,2) = (xM , yM , zM ) + α1,2(nx, ny, nz). If D = 0, then M1 ≡M2; for D < 0 the line does not intersect the sphere. The Cartesian coordinates of intersections need to be converted to the spherical, where ϕ1,2 = arcsin ( z1,2/ √ x2 1,2 + y2 1,2 + z2 1,2 ) , λ1,2 = atan2(y1,2, x1,2). Practical computation. The meridian plane %1 containing m(λc) passes through the points P5 = [ϕ5, λc], P6 = [ϕ6λc], and P7 = [ϕ7, λc], where ϕ5 = −π 2 + ε, ϕ6 = 0, ϕ7 = π 2 − ε. Analogously, the parallel plane %2 containing p(ϕc, λ) is given by the points P8 = [ϕc, λ8], P9 = [ϕc, λ9], and P10 = [ϕc, λ10], where λ8 = −π + ε, λ9 = 0, λ10 = π − ε, and the meridian plane %3 containing m′(ϕ′, λ′0) centered at λ′0 is given by three points P11 = [ϕ′11, λ ′ 0], P12 = [ϕ′12, λ ′ 0], and P13 = [ϕ′13, λ ′ 0], where ϕ′11 = −π 2 + ε, ϕ′12 = 0, ϕ′13 = π 2 − ε. 1. The inverse transformation Initially, the spherical coordinates [ϕ′i, λ′i] of points P11, P12, P13 related to the pole K = [ϕk, λk] are transformed to the normal aspect using Eqs. 3.4, 3.5, where i = 11, 12, 13. 2. The transformation to the geocentric coordinates The spherical coordinates [ϕi, λi] of points P5 − P13 are transformed to the Cartesian geocentric coordinates xi, yi, zi, where xi = cosϕi cosλi, yi = cosϕi sinλi, zi = sinϕi, i = 5, ..., 13. 3. Find the analytic representation of planes The meridian plane of m′(λ′0) is given by the vectors u1, v1, where u1 = P12−P11, v1 = P13−P11. The vectors u2, v2 of the second meridian plane m(λc) are u2 = P6−P5, v2 = P7−P5, and the vectors u3, v3 of the parallel plane p(ϕc) are u3 = P9−P8, v3 = P10−P8. Find the vectors n1, n2, n3 perpendicular to u1, v1, u2, v2, and u3, v3 using the cross products n1 = u1 × v1, n2 = u2 × v2, n3 = u3 × v3. 4. Compute the line of intersection For the meridian planes m′(λ′0), and m(λc), set: Q1 ≡ P11, Q2 ≡ P5; for the meridian m′(λ′0) and parallel p(ϕc) planes, set: Q1 ≡ P11, Q2 ≡ P8, compute Q0 = (P11 + P12 + P13)/3. Evaluate the line of intersection, compute the cross products n1,2 = n1 × n2, n1,3 = n1×n3, where the direction n of both lines of intersection is n ≡ n1,2, or n ≡ n1,3. 5. Intersection of the line and sphere Let S2 be the sphere of the radius R = 1 centered at C = [0, 0, 0]. Depending on the sign of the discriminantD, evaluate the intersectionsM1 = [x1, y1, z1],M2 = [x2, y2, z2] of the Geoinformatics FCE CTU 17(2), 2018 39 T. Bayer: Plotting the map projection graticule with discontinuities -18 0.0 -15 0.0 -12 0.0 -90 .0 -60 .0 -30.0 0.030.0 60.0 90 .0 12 0.0 15 0.0 18 0.0 15.015.0 -16 5.0 -16 5.0 -90.0 -90.0-90.0 -75.0 -75.0 -75.0 -60.0 -60.0 -60.0 -45.0 -45.0 -45.0 -30.0 -30.0 -30.0 -15.0 -15 .0 -15.0 0.0 0.0 0.0 15.0 15 .0 15.0 30.0 30 .0 30.0 45.0 45 .0 45.0 60.0 60 .0 60.0 75.0 75 .0 75.0 90.090 .0 90.0 -18 0.0 -15 0.0 -12 0.0 -90 .0 -60 .0 -30.0 0.030.0 60.0 90 .0 12 0.0 15 0.0 18 0.0 15.015.0 -16 5.0 -16 5.0 -90.0 -75.0 -60.0 -45.0 -30 .0 -15 .0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 a) b) Figure 4.1: The Hassler projection in the oblique aspect, K = [50◦, 15◦]; the singularities involved a) and treated b). planes %2 and %3. The planes %1, %3 have two intersections denoted as M3 = [x3, y3, z3] and M4 = [x4, y4, z4]. 6. Transformation to spherical coordinates Finally, the Cartesian coordinates of the intersections M1, M2, M3 will be converted to the spherical. 4. Graticule reconstruction using the combined sampling Let Ω = Ωϕ × Ωλ, where Ωϕ = [ϕ,ϕ], and Ωλ = [λ, λ], be the subdomain, inside which the polygonal approximation of the meridians/parallels is constructed. The points pi = (F (ϕi, λc), G(ϕi, λc)), 1 ≤ i ≤ nϕ, ϕi ∈ Ωϕ approximate the meridian m(λc), 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 (pi−1, pi, pi+1). 4.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 qi = [ϕi, λi] is “good” if no singularity at F (qi) and G(qi) 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. Geoinformatics FCE CTU 17(2), 2018 40 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 1 Combined sampling with the singularities, the stack implementation. 1: function Graticule(Lϕ,Lλ, Ωϕ, Ωλ,∆ϕ,∆λ,P, s, d, d, ε, α) 2: S = ∅, s = 0 3: S ← Ωj = Ωϕ ×Ωλ 4: while S 6= ∅ do 5: Ωj = Ωϕ,j ×Ωλ,j = S.pop() 6: Ωj,1 = Ωϕ,j,1 ×Ωλ,j,1 = Ωj , Ωj,2 = Ωϕ,j,2 ×Ωλ,j,2 = Ωj 7: k = 0 8: try 9: Lϕ,j = ∅, Lλ,j = ∅ 10: Meridians(Lλ,j , Ωϕ,j , Ωλ,j ,∆λ,P, d, d, ε, α) 11: Parallels(Lϕ,j , Ωϕ,j , Ωλ,j ,∆ϕ,P, d, d, ε, α) 12: Lϕ ← Lϕ,j , Lλ ← Lλ,j 13: catch (SingularityLatError e) 14: if (s < s) then 15: processInt(Ωϕ,j , e.ϕi, Ωϕ,j,1, Ωϕ,j,2, k, s, ε) 16: catch (SingularityLonError e) 17: if (s < s) then 18: processInt(Ωλ,j , e.λi, Ωλ,j,1, Ωλ,j,2, k, s, ε) 19: catch (SingularityError e) 20: if (s > s) then 21: Lϕ = ∅,Lλ = ∅ 22: return 23: else 24: if (k > 0) then 25: S ← Ωj,1 26: else if (k > 1) then 27: S ← Ωj,2 Initially, Ωg ≡ Ω is set, so Ωg ϕ ≡ Ωϕ, Ωg ϕ ≡ Ωϕ. All adaptively sampled meridian/parallel points qi 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: Ωϕ,j = 〈Ωϕ,j,1, Ωϕ,j,2〉 , Ωϕ,j,1 = [ϕ j , c− ε], Ωϕ,j,2 = [c+ ε, ϕj ], Ωλ,j = 〈Ωλ,j,1, Ωλ,j,2〉 , Ωλ,j,1 = [λj , c− ε], Ωλ,j,2 = [c+ ε, λj ]. Otherwise, the polygonal approximation Lλ,j of the meridian m(λ) or Lϕ,j of the parallel p(ϕ) is constructed. All disjoint polygonal approximations Lj 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. Geoinformatics FCE CTU 17(2), 2018 41 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 2 The construction of meridians. 1: function Meridians(Lλ,j , Ωϕ,j , Ωλ,j ,∆λ,P, d, d, ε, α) 2: λj,start = ∆λ[1 + floor( λj , ∆λ)] 3: λj,end = ∆λ[ceil( λj ∆λ)− 1] 4: meridianPart(Lλ,j , Ωϕ,j , λj,start,P, d, d, ε, α) 5: for (λ = λj,start;λ ≤ λj,end;λ = λ+ ∆λ) do 6: meridianPart(Lλ,j , Ωϕ,j , λ,P, d, d, ε, α) 7: if (λk ∈ Ωλ,j) then 8: meridianPart(Lλ,j , Ωϕ,j , λ− ε,P, d, d, ε, α) 9: meridianPart(Lλ,j , Ωϕ,j , λ+ ε,P, d, d, ε, α) 10: λk,o = (λk > 0?λk − π : λk + π) 11: meridianPart(Lλ,j , Ωϕ,j , λk,o − ε,P, d, d, ε, α) 12: meridianPart(Lλ,j , Ωϕ,j , λk,o + ε,P, d, d, ε, α) 13: meridianPart(Lλ,j , Ωϕ,j , λj,end,P, d, d, ε, α) 2. 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 Ωϕ,j = [ϕ j , ϕj ], Ωλ,j = [λj , λj ]. (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Ωj,2 = Ωϕ,j,2×Ωλ,j,2, as Ωj,1 = Ωj , Ωj,2 = Ωj , where Ωϕ,j,1 = [ϕ j,1, ϕj,1], Ωλ,j,1 = [λj,1, λj,1], and Ωϕ,j,2 = [ϕ j,2, ϕj,2], Ωλ,j,2 = [λj,2, λj,2] 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 toLλ: Lλ ← Lλ,j , and Lϕ,j toLϕ: 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). (f) Resolve the singularities by splitting Ωϕ,j , Ωλ,j Geoinformatics FCE CTU 17(2), 2018 42 T. Bayer: Plotting the map projection graticule with discontinuities Depending on the singularity c direction, do the following steps: i. If c = ϕi, then Ωϕ,j needs to be split. Call the function processInt() with parameters: Ωϕ,j , c, Ωϕ,j,1, Ωϕ,j,2, k, S, ε. The subintervals Ωϕ,j,1, Ωϕ,j,2 are passing by reference. ii. If c = λi, then Ωλ,j needs to be split. Call the function processInt() with parameters: Ωλ,j , c, Ωλ,j,1, Ωλ,j,2, k, S, ε. The subintervals Ωλ,j,1, Ωλ,j,2 are passing by reference. 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 fail- ures. The singularity value c is stored in the thrown exceptions SingularityLatError, SingularityLonError, which are derived from the base class SingularityError. The split- ting procedure is realized in accordance with the singularity direction. Generate meridians. The procedure creates all meridians inside the given interval Ωj = Ωϕ,j × Ωλ,j , where Ωϕ,j = [ϕ j , ϕj ], Ωλ,j = [λj , λj ]. 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 λj,start = ∆λ[1 + floor( λj , ∆λ)], λj,end = ∆λ[ceil( λj∆λ)− 1]. 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 Ωj = Ωϕ,j × Ωλ,j , where Ωϕ,j = [ϕ j , ϕj ], Ωλ,j = [λj , λj ]. 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 meridianm(λ) by the refinement α and the recursion depth d over the interval Ωϕ,j = [ϕ j , ϕj ]. Geoinformatics FCE CTU 17(2), 2018 43 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 3 The construction of parallels. 1: function Parallels(Lλ,j , Ωϕ,j , Ωλ,j ,∆ϕ,P, d, d, ε, α) 2: ϕj,start = ∆ϕ[1 + floor( ϕ j , ∆ϕ)] 3: ϕj,end = ∆ϕ[ceil( ϕj ∆ϕ)− 1] 4: parallelPart(Lϕ,j , Ωϕ,j , ϕj,start,P, d, d, ε, α) 5: for (ϕ = ϕj,start;ϕ ≤ ϕj,end;ϕ = ϕ+ ∆ϕ) do 6: parallelPart(Lϕ,j , Ωϕ,j , ϕ,P, d, d, ε, α) 7: parallelPart(Lϕ,j , Ωϕ,j , ϕj,end,P, d, d, ε, α) Suppose m′(λ′0) to be the central meridian. The polygonal approximation consists of two steps: 1. The partition of m(λ) to the fragments. Depending on the position of the pole K = [ϕk, λk] and the central meridian m′(λ′0), the original interval Ωϕ,j of the meridian m(λ) needs to be split along the intersection point M3 = [ϕ3, λ]. If λ′0 = 0 the meridians m′(0) and m(λ) have just two intersections M3 = [π/2, ·] ≡ A, M4 = [−π/2, ·] ≡ B at poles A,B (this case is not important for the algorithm). Otherwise, the solution M3 = [ϕ3, λ3] found from the intersection of two planes (see Sec. 3.2) holds ϕ3 ∈ [ϕ j + ε, ϕj − ε] ∧ (|λ3 − λ| < ε) . If an intersection M3 is found, the splitting procedure Ωϕ,j,1 = [ϕ j , ϕ3 − ε], Ωϕ,j,2 = [ϕ3 + ε, ϕj ], is undertaken; see Fig. 3.4 and Alg. 4. 2. 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 meridianm′(λ′0) in at most two points. The original interval Ωλ,j needs to be split along the intersection pointsM1 = [ϕ1, λ1], M2 = [ϕ2, λ2]; the solution is found from the intersection of two planes (see Sec. 3.2). In general, the following cases are solved: • If M1 ∈ Ωλ,j , than λ1 ∈ [λj + ε, λj − ε] and Ωλ,j is split to Ωλ,j,1 = [λj , λ1 − ε], Ωλ,j,2 = [λ1 + ε, λj ]. • If M2 ∈ Ωλ,j , than λ2 ∈ [λj + ε, λj − ε] and Ωλ,j is split to Ωλ,j,1 = [λj , λ2 − ε], Ωλ,j,2 = [λ2 + ε, λj ]. • If M1 ∈ Ωλ,j ∧M1 ∈ Ωλ,j , than Ωλ,j is split to Ωλ,j,1 = [λj , λ1 − ε], Ωλ,j,2 = [λ1 + ε, λ2 − ε], Ωλ,j,3 = [λ2 + ε, λj ]. For practical computation, the coordinates are sorted so that λ1 ≤ λ2; see Alg. 5. Geoinformatics FCE CTU 17(2), 2018 44 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 4 The construction of a meridian from the fragments. 1: function meridianPart(Lλ,j,k, Ωϕ,j , λ,P, d, d, ε, α) 2: 〈 M1 3 ,M 2 3 〉 = m(λ)×m′(ϕ′, λ′0) 3: Lλ,j,1 = ∅, Lλ,j,2 = ∅ 4: if ( ϕ1 3 > ϕ j + ε ) ∧ ( ϕ1 3 < ϕj − ε ) ∧ (∣∣λ1 3 − λ ∣∣ < ε ) then 5: ϕj,1 = ϕ1 3 − ε, ϕj,2 = ϕ1 3 + ε 6: csMerInit(P, Lλ,j,1, ϕj , ϕj,1, λ, d, d, ε, α) 7: Lλ,j ← Lλ,j1 8: csMerInit(P, Lλ,j,2, ϕj,2, ϕj , λ, d, d, ε, α) 9: Lλ,j ← Lλ,j2 10: else if ( ϕ2 3 > ϕ j + ε ) ∧ ( ϕ2 3 < ϕj − ε ) ∧ (∣∣λ2 3 − λ ∣∣ < ε ) then 11: ϕj,1 = ϕ2 3 − ε, ϕj,2 = ϕ2 3 + ε 12: csMerInit(P, Lλ,j,1, ϕj , ϕj,1, λ, d, d, ε, α) 13: Lλ,j ← Lλ,j1 14: csMerInit(P, Lλ,j,2, ϕj,2, ϕj , λ, d, d, ε, α) 15: Lλ,j ← Lλ,j2 4.2. 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. 1. The initial phase Let Lλ,j,k = ∅ be the empty set. Compute xa = F (a, λ), ya = G(a, λ) and xb = F (b, λ), yb = G(b, λ), where pa = [xa, ya] and pb = [xb, yb]. If a singularity in a or b is detected, throw an exception. Add the initial vertex to Lλ,j,k: Lλ,j,k ← pa, where pa = (xa, ya). Set the recursion depth d = 1. The procedure is summarized in Alg. 6. 2. The recursive procedure Enter the recursive procedure and do the following substeps: (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 ϕ1 = a+ 1 2r1(b− a), ϕ2 = a+ r2(b− a), ϕ3 = a+ 3 2r3(b− a), into the approximate quarters Ωϕ,j,k,1 = [a, ϕ1], Ωϕ,j,k,2 = [ϕ1, ϕ2] Ωϕ,j,k,3 = [ϕ2, ϕ3], Ωϕ,j,k,4 = [ϕ3, b], where r1, r2, r3 are the random numbers inside the interval [0.45, 0.55]. Geoinformatics FCE CTU 17(2), 2018 45 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 5 The construction of a parallel from the fragments. 1: function parallelPart(Lϕ,j,k, Ωλ,j , ϕ,P, d, d, ε, α) 2: 〈M1,M2〉 = p(ϕ)×m′(ϕ′, λ′0) 3: Lϕ,j,1 = ∅, Lϕ,j,2 = ∅, Lϕ,j,3 = ∅ 4: if (λ1 > λj + ε) ∧ (λ1 < λj − ε) ∧ (λ2 > λj + ε) ∧ (λ2 < λj − ε) then 5: λ1,2 = min(λ1, λ2), λ1,2 = max(λ1, λ2) 6: λj,1 = λ1,2 − ε, λj,2 = λ1,2 + ε, λj,2 = λ1,2 − ε, λj,2 = λ1,2 + ε 7: csParInit(P, Lϕ,j,1, λj , λj,1, ϕ, d, d, ε, α) 8: Lϕ,j ← Lϕ,j1 9: csParInit(P, Lϕ,j,2, λj,2, λj,2, ϕ, d, d, ε, α) 10: Lϕ,j ← Lϕ,j2 11: csParInit(P, Lϕ,j,3, λj,3, λj , ϕ, d, d, ε, α) 12: Lϕ,j ← Lϕ,j3 13: else if (λ1 > λj + ε) ∧ (λ1 < λj − ε) then 14: λj,1 = λ1 − ε, λj,2 = λ1 + ε 15: csParInit(P, Lϕ,j,1, λj , λj,1, ϕ, d, d, ε, α) 16: Lϕ,j ← Lϕ,j1 17: csParInit(P, Lϕ,j,2, λj2 , λj , ϕ, d, d, ε, α) 18: Lϕ,j ← Lϕ,j2 19: else if (λ2 > λj + ε) ∧ (λ2 < λj − ε) then 20: λj,1 = λ2 − ε, λj,2 = λ2 + ε 21: csParInit(P, Lϕ,j,1, λj , λj,1, ϕ, d, d, ε, α) 22: Lϕ,j ← Lϕ,j1 23: csParInit(P, Lϕ,j,2, λj2 , λj , ϕ, d, d, ε, α) 24: Lϕ,j ← Lϕ,j2 (c) If a singularity in [ϕ1, λ], [ϕ2, λ], or [ϕ3, λ] occurs, throw a new exception according to the discontinuity direction. (d) Otherwise, evaluate the function values x1 = F (ϕ1, λ), y1 = G(ϕ1, λ), x2 = F (ϕ2, λ), y2 = G(ϕ2, λ), x3 = F (ϕ3, λ), y3 = G(ϕ3, λ), at new vertices p1, p2, p3. (e) 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 α1 = α(pa, p1, p2), α2 = α(p1, p2, p3), α3 = α(p2, p3, pb), and the recursive depth d. When the meridian is not sufficiently smooth, or d ≤ d, it needs to be refined. Hence, the recursive subdivision is necessary. (f) If α1 > α, call the recursive procedure with the increased depth d = d+ 1 for the interval [a, ϕ1]. (g) Add the new point p1 to the polynomial approximation of m̃(λ) : Lλ,j,k ← p1. (h) If α1 > α∨α2 > α, call the recursive procedure with the increased depth d = d+1 for the interval [ϕ1, ϕ2]. (i) Add the new point p2 to the polynomial approximation of m̃(λ) : Lλ,j,k ← p2. Geoinformatics FCE CTU 17(2), 2018 46 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 6 Combined sampling of the meridian, the initial phase. 1: function csMerInit(P, Lλ,j,k, a, b, λ, d, d, d, ε, α) 2: Lλ,j,k = ∅ 3: if discontinuity in (a, λ) in ϕ direction then 4: throw LatSingularityException (a) 5: if discontinuity in (a, λ) in λ direction then 6: throw LonSingularityException (λ) 7: if discontinuity in (b, λ) in ϕ direction then 8: throw LatSingularityException (b) 9: if discontinuity in (b, λ) in λ direction then 10: throw LonSingularityException (λ) 11: xa = F (a, λ), ya = G(a, λ) 12: Lλ,j,k ← Point(xa, ya) 13: csMer(P, Lλ,j,k, a, b, xa, ya, xb, yb, d, d, d, ε, α) 14: xb = F (b, λ), yb = G(b, λ) 15: Lλ,j,k ← Point(xb, yb) (j) If α2 > α∨α3 > α, call the recursive procedure with the increased depth d = d+1 for the interval [ϕ2, x3]. (k) Add the new point p3 to the polynomial approximation of m̃(λ) : Lλ,j,k ← p3. (l) If α3 > α, call the recursive procedure for the interval [ϕ3, b]. 3. Add the last point Add the last point pb to the polynomial approximation of m̃(λ) : Lλ,j,k ← pb 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 α1 > α, α1 > α ∨ α2 > α, α2 > α ∨ α3 > α, and α3 > α, are replaced with the new α1 > α ∨ α2 > α ∨ α3 > α, (4.1) 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 ac- count angles between the newly created and adjacent segments from the previous iteration. Geoinformatics FCE CTU 17(2), 2018 47 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 7 Combined sampling of the meridian points, the recursive phase, method S1. 1: function csMer(P, Lλ,j,k, a, b, xa, ya, xb, yb, d, d, d, ε, α)) 2: if (d > d) ∨ (b− a < ε) then 3: return 4: r1 = rand(0.45, 0.55), r2 = rand(0.45, 0.55), r3 = rand(0.45, 0.55) 5: ϕ1 = a+ 1 2r1(b− a), ϕ2 = a+ r2(b− a), ϕ3 = a+ 3 2r3(b− a) 6: if discontinuity in (ϕi, λ)i in ϕ direction then 7: throw LatSingularityException (ϕi), i = 1, 2, 3 8: if discontinuity in (ϕi, λ)i in λ direction then 9: throw LonSingularityException (λ), i = 1, 2, 3 10: xi = F (ϕi, λ), yi = G(ϕi, λ), i = 1, 2, 3 11: pa = Point(xa, ya), pb = Point(xb, yb), pi = Point(xi, yi), i = 1, 2, 3 12: α1 = α(pa, p1, p2), α2 = α(p1, p2, p3), α3 = α(p2, p3, pb) 13: if (α1 > α) ∨ (d <= d) then 14: csMer(P, Lλ,j,k, a, ϕ1, xa, ya, x1, y1, d+ 1, d, d, α); 15: Lλ,j,k ← Point(x1, y1) 16: if (α1 > α) ∨ (α2 > α) ∨ (d <= d) then 17: csMer(P, Lλ,j,k, ϕ1, ϕ2, x1, y1, x2, y2, d+ 1, d, d, α); 18: Lλ,j,k ← Point(x2, y2) 19: if (α2 > α) ∨ (α3 > α) ∨ (d <= d) then 20: csMer(P, Lλ,j,k, ϕ2, ϕ3, x2, y2, x3, y3, d+ 1, d, d, α); 21: Lλ,j,k ← Point(x3, y3) 22: if (α3 > α) ∨ (d <= d) then 23: csMer(P, Lλ,j,k, ϕ3, b, x3, y3, xb, yb, d+ 1, d, d, α); For a given Ωϕ,j,k = [ak, bk], its predecessor Ωϕ,j,k−1 = [ak−1, bk−1] and successor Ωϕ,j,k+1 = [ak+1, bk+1] are Ωϕ,j,k−1 = pred(Ωϕ,j,k), Ωϕ,j,k+1 = succ(Ωϕ,j,k), where bk−1 = ak, ak+1 = bk. Analogously, for the polygonal approximation of Ωϕ,j,k bounds are pa,k−1 = pred(pa,k), pb,k+1 = succ(pb,k), where pa,k−1 = [xa,k−1, ya,k−1], pa,k = [xa,k, ya,k] ≡ pb,k−1, pb,k = [xb,k, yb,k] ≡ pa,k+1, pb,k+1 = [xb,k+1, yb,k+1], and xa,k−1 = F (ak−1, λ), xa,k = F (ak, λ), xb,k = F (bk, λ), xb,k+1 = F (bk+1, λ); see Fig. 4.2 The approximate quarters of Ωϕ,j,k denoted as Ωϕ,j,k,l, l = 1, ..., 4, have the form of Ωϕ,j,k,1 = [ak, ϕ1,k], Ωϕ,j,k,2 = [ϕ1,k, ϕ2,k] Ωϕ,j,k,3 = [ϕ2,k, ϕ3,k], Ωϕ,j,k,4 = [ϕ3,k, bk], The idea of modification is straightforward, an information about previous and subsequent quarters of Ωϕ,j,k,l is stored Ωϕ,j,k,l−1 = pred(Ωϕ,j,k,l), Ωϕ,j,k,l+1 = succ(Ωϕ,j,k,l). While the first quarter lacks its predecessor, the last quarter lacks the successor. For a current recursive depth d, the segment (pa,k, pb,k) of the polygonal approximation is split to the quarters, where p1,k, p2,k, p3,k are the newly created vertices at the recursive depth d+ 1 representing the “projected bounds” of newly created intervals Ωϕ,j,k,1,...,Ωϕ,j,k,4. Unlike S1, S2, the angles αi are measured not only between the newly created segments Geoinformatics FCE CTU 17(2), 2018 48 T. Bayer: Plotting the map projection graticule with discontinuities p a,1 p b,1 m(l) p a,1 p b,1 =p a,2 p 1,1 p 2,1 p 3,1 p 1,1 p 2,1 p 3,1 p 1,2 p 2,2 p 3,2 p 1,3 p 2,3 p 3,3 m(l) m(l) m(l) p 3,1 p 1,3 p a,1 p b,4 p 3,2 p 1,4 p a,1 p b,1 =p a,2 p b,2 =p a,3 p b,3 =p a,4 p b,1 =p a,2 p b,2 =p a,3 p b,3 =p a,4 p b,4 p b,4 p b,2 =p a,3 p b,3 =p a,4 p 1,2 Figure 4.2: Illustration of the combined sampling method S3 for the meridian m(λ). α1 = α(pa,k, p1,k, p2,k), α2 = α(p1,k, p2,k, p3,k), α3 = α(p2,k, p3,k, pb,k), in the recursive depth d+ 1, but also between the adjacent segments α0 = α(p3,k−1, pa,k, p1,k), α4 = α(p3,k, pb,k, p1,k+1), where p3,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 ϕ3,k−1 = ak + 1 2r0(ak−1 − ak), p1,k+1 is the projected first-quarter Ωϕ,j,k+1,1 of the next interval Ωϕ,j,k+1 with the latitude ϕ1,k+1 = bk + 1 2r4(bk+1 − bk), and r0, r4 are the random numbers generated according to the principles mentioned above. The modified condition for the recursive call has the form of α0 > α ∨ α1 > α ∨ α2 > α ∨ α3 > α ∨ α4 > α. (4.2) 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: 1. If d > d or bk − ak < ε, stop the recursive procedure and go to Step 3. 2. For a given Ωϕ,j,k = [ak, bk], its predecessor Ωϕ,j,k−1 = [ak−1, bk−1], and successor Ωϕ,j,k+1 = [ak+1, bk+1], the interval is split by three points, the positions are close to the quarters of the interval ϕ1,k = ak + 1 2r1(bk − ak), ϕ2,k = ak + r2(bk − ak), ϕ3,k = ak + 3 2r3(bk − ak). Geoinformatics FCE CTU 17(2), 2018 49 T. Bayer: Plotting the map projection graticule with discontinuities Algorithm 8 Combined sampling of the meridian points, the recursive phase, method S3. 1: function csMer(P, Lλ,j,k, ak, bk, ak−1, bk+1, xa, ya, xb, yb, d, d, d, ε, α)) 2: if (d > d) ∨ (bk − ak < ε) then 3: return 4: ri = rand(0.45, 0.55), i = 0, ..., 4 5: ϕ1,k = ak + 1 2r1(bk − ak), ϕ2,k = a+ r2(bk − a), ϕ3,k = a+ 3 2r3(bk − ak) 6: ϕ3,k−1 = ak + 1 2r1(bk−1 − ak), ϕ1,k+1 = bk + r4(bk+1 − bk) 7: if discontinuity in (ϕi,k, λ) in ϕ direction then 8: throw LatSingularityException (ϕi,k), i = 1, 2, 3 9: if discontinuity in (ϕi,k, λ) in λ direction then 10: throw LonSingularityException (λ), i = 1, 2, 3 11: x3,k−1 = F (ϕ3,k−1, λ), y3,k−1 = G(ϕ3,k−1, λ) 12: xi,k = F (ϕi,k, λ), yi,k = G(ϕi,k, λ), i = 1, 2, 3 13: x1,k+1 = F (ϕ1,k+1, λ), y1,k+1 = G(ϕ1,k+1, λ) 14: pa,k = Point(xa,k, ya,k), pb,k = Point(xb,k, yb,k) 15: pi,k = Point(xi,k, yi,k), i = 1, 2, 3 16: p3,k−1 = Point(x3,k−1, y3,k−1), p1,k+1 = Point(x1,k+1, y1,k+1) 17: α1 = α(pa,k, p1,k, p2,k), α2 = α(p1,k, p2,k, p3,k), α3 = α(p2,k, p3,k, pb,k) 18: α0 = α(p3,k−1, pa,k, p1,k), α4 = α(p3,k, pb,k, p1,k+1) 19: if (α0 > α) ∨ (α1 > α) ∨ (α2 > α) ∨ (α3 > α) ∨ (α4 > α) ∨ (d <= d)) then 20: csMer(P, Lλ,j,k, ak, ϕ1,k, ϕ3,k−1, ϕ2,k, xa,k, ya,k, x1,k, y1,k, d+ 1, d, d, α) 21: Lλ,j,k ← Point(x1,k, y1,k) 22: if (α0 > α) ∨ (α1 > α) ∨ (α2 > α) ∨ (α3 > α) ∨ (α4 > α) ∨ (d <= d) then 23: csMer(P, Lλ,j,k, ϕ1,k, ϕ2,k, ϕa,k, ϕ3,k, x1,k, y1,k, x2,k, y2,k, d+ 1, d, d, α) 24: Lλ,j,k ← Point(x2,k, y2,k) 25: if (α0 > α) ∨ (α1 > α) ∨ (α2 > α) ∨ (α3 > α) ∨ (α4 > α) ∨ (d <= d) then 26: csMer(P, Lλ,j,k, ϕ2,k, ϕ3,k, ϕ1,k, bk, x2,k, y2,k, x3,k, y3,k, d+ 1, d, d, α) 27: Lλ,j,k ← Point(x3,k, y3,k) 28: if (α0 > α) ∨ (α1 > α) ∨ (α2 > α) ∨ (α3 > α) ∨ (α4 > α) ∨ (d <= d) then 29: csMer(P, Lλ,j,k, ϕ3,k, bk, ϕ2,k, ϕ1,k+1, x3,k, y3,k, xb,k, yb,k, d+ 1, d, d, α) Subsequently, the third quarter ϕ3,k−1 of the previous interval Ωϕ,j,k−1 representing the lower bound of Ωϕ,j,k−1,4 and the first quarter ϕ1,k+1 of the next intervals Ωϕ,j,k+1 representing the upper bound of Ωϕ,j,k+1,1, are evaluated. (a) If a singularity in [ϕ3,k−1, λ], [ϕ1,k, λ], [ϕ2,k, λ], [ϕ3,k, λ], or [ϕ1,k+1, λ] occurs, throw a new exception according to the discontinuity direction. (b) Otherwise, evaluate the function values x3,k−1 = F (ϕ3,k−1, λ), y3,k−1 = G(ϕ3,k−1, λ), x1,k = F (ϕ1,k, λ), y1,k = G(ϕ1,k, λ), x2,k = F (ϕ2,k, λ), y2,k = G(ϕ2, λ), x3,k = F (ϕ3,k, λ), y3,k = G(ϕ3,k, λ), x1,k+1 = F (ϕ1,k+1, λ), y1,k+1 = G(ϕ1,k+1, λ), at new vertices p1,k, p2,k, p3,k and the supplementary vertices p3,k−1, p1,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 α0 = α(p3,k−1, pa,k, p1,k), α1 = α(pa,k, p1,k, p2,k), α2 = α(p1,k, p2,k, p3,k), α3 = Geoinformatics FCE CTU 17(2), 2018 50 T. Bayer: Plotting the map projection graticule with discontinuities Figure 4.3: The discontinuity treatment without its classification using the proposed algorithm in polyconic projection; the doubled equator is visible under magnification. α(p2,k, p3,k, pb,k), α4 = α(p3,k, pb,k, p1,k+1), 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 d = d + 1 for the interval Ωϕ,j,k,1 = [ak, ϕ1,k], its predecessor Ωϕ,j,k−1,4 = [ϕ3,k−1, ak], and successor Ωϕ,j,k,2 = [ϕ1,k, ϕ2,k]. (e) Add the new point p1,k to the polynomial approximation of m̃(λ) : Lλ,j,k ← p1,k. (f) If the condition given by Eq. 4.2 is held, call the recursive procedure with the increased depth d = d + 1 for the interval Ωϕ,j,k,2 = [ϕ1,k, ϕ2,k], its predecessor Ωϕ,j,k,1 = [ak, ϕ1,k], and successor Ωϕ,j,k,3 = [ϕ2,k, ϕ3,k]. (g) Add the new point p2,k to the polynomial approximation of m̃(λ) : Lλ,j,k ← p2,k. (h) If the condition given by Eq. 4.2 is held, call the recursive procedure with the increased depth d = d + 1 for the interval Ωϕ,j,k,3 = [ϕ2,k, ϕ3,k], its predecessor Ωϕ,j,k,2 = [ϕ1,k, ϕ2,k], and successor Ωϕ,j,k,4 = [ϕ3,k, bk]. (i) Add the new point p3,k to the polynomial approximation of m̃(λ) : Lλ,j,k ← p3,k. (j) If the condition given by Eq. 4.2 is held, call the recursive procedure for the last interval Ωϕ,j,k,4 = [ϕ3,k, bk], its predecessor Ωϕ,j,k,3 = [ϕ2,k, ϕ3,k], and successor Ωϕ,j,k+1,1 = [bk, ϕ1,l+1]. 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 Geoinformatics FCE CTU 17(2), 2018 51 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 λ1 = a+ 1 2r1(b− a), λ2 = a+ r2(b− a), λ3 = a+ 3 2r3(b− a), split the interval Ωλ,j,k to the approximate quarters Ωλ,j,k,1 = [a, λ1], Ωλ,j,k,2 = [λ1, λ2] Ωλ,j,k,3 = [λ2, λ3], Ωλ,j,k,4 = [λ3, b]. 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 disconti- nuity 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. 5. Experiments and results The principles mentioned above will be illustrated on the several tests. While the first analy- sis 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 coor- dinate functions are involved). Finally, the last experiment illustrates the splitting procedure efficiency depending on the threshold ε. The results are summarized in Tabs. 1-5. 5.1. 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 nmer and parallel points npar. 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 satis- fying 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 sam- pling 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 nmer, npar is analogous. In Geoinformatics FCE CTU 17(2), 2018 52 T. Bayer: Plotting the map projection graticule with discontinuities Table 1: Combined sampling of the graticule for α = 1◦, 2◦, 5◦, methods S1, S2, S3 are compared. Projection # nmer npar αm αp α̃m α̃p Equidistant cylindrical S1 195/195/195 190/190/190 0/0/0 0/0/0 0/0/0 0/0/0 S2 195/195/195 190/190/190 0/0/0 0/0/0 0/0/0 0/0/0 S3 195/195/195 190/190/190 0/0/0 0/0/0 0/0/0 0/0/0 Equidistant conic S1 195/195/195 9961/6973/2464 0/0/0 1.02/2.27/4.98 0/0/0 0.53/0.75/2.16 S2 195/195/195 10066/7678/2470 0/0/0 1.00/2.13/4.06 0/0/0 0.52/0.69/2.16 S3 195/195/195 10090/8734/2470 0/0/0 0.99/1.99/4.41 0/0/0 0.52/0.60/2.16 Equidistant azimuthal S1 195/195/195 12352/9091/2494 0/0/0 1.08/2.04/4.87 0/0/0 0.56/0.76/2.78 S2 195/195/195 13186/9238/2494 0/0/0 1.09/2.06/4.89 0/0/0 0.52/0.74/2.78 S3 195/195/195 13786/9706/2506 0/0/0 1.00/1.90/4.78 0/0/0 0.50/0.71/2.77 Bonne S1 14250/7206/2901 6610/3613/1576 1.21/2.27/5.68 1.05/2.06/4.98 0.45/0.88/2.20 0.46/0.84/1.98 S2 15171/7599/3099 6766/3778/1630 1.05/2.52/5.00 1.00/2.00/4.82 0.42/0.84/2.07 0.45/0.81/1.91 S3 16839/8367/3483 6958/4330/1654 1.00/1.99/4.98 1.00/1.99/4.96 0.38/0.76/1.84 0.44/0.71/1.82 Hassler S1 16452/7170/3342 8839/5779/1912 1.14/2.13/5.77 1.22/2.25/4.91 0.41/0.96/2.09 0.51/0.78/2.39 S2 17391/7479/3519 9370/5974/1966 1.07/2.17/5.20 1.03/2.11/4.99 0.40/0.92/1.99 0.48/0.75/2.32 S3 19899/8379/4191 9826/6670/2002 1.00/2.00/4.96 1.00/1.99/4.76 0.35/0.82/1.66 0.46/0.67/2.28 Nicolosi S1 14018/6503/3227 5519/2495/1124 1.14/2.29/5.64 1.15/2.16/5.50 0.43/0.94/1.92 0.42/0.95/2.12 S2 14468/6944/3320 5804/2564/1208 1.00/2.00/4.96 1.01/2.00/4.94 0.42/0.88/1.86 0.40/0.92/1.96 S3 15176/8156/3488 6224/2912/1352 1.00/2.00/4.93 1.00/1.99/4.89 0.40/0.74/1.77 0.37/0.81/1.76 Littrow S1 6513/3279/1470 12011/7052/2507 1.10/2.15/5.43 1.14/2.27/5.17 0.50/1.00/2.30 0.51/0.87/2.48 S2 6963/3615/1587 12242/7334/2606 1.06/2.22/4.90 1.05/1.99/5.27 0.47/0.91/2.12 0.50/0.84/2.38 S3 8121/4374/2151 13796/7931/2903 1.00/2.00/4.86 1.00/1.98/4.98 0.40/0.75/1.55 0.45/0.77/2.13 Adams I. S1 7083/3669/1788 4430/2285/887 1.07/2.00/4.93 1.06/2.10/6.16 0.44/0.86/1.82 0.44/0.83/2.04 S2 7449/4017/2097 4778/2534/1046 1.07/2.08/4.96 1.02/2.01/4.86 0.42/0.78/1.53 0.41/0.74/1.74 S3 8085/4689/2148 4562/2618/1058 1.00/2.00/4.89 0.98/1.99/4.44 0.38/0.67/1.49 0.40/0.72/1.72 Table 2: Uniform sampling of the graticule with the determined steps δϕ, δλ compared to S3; the amounts of sampled points nmer, npar are preserved. Projection α δϕ δλ nmer ∆nmer npar ∆npar αm ∆αm αp ∆αp α̃m ∆α̃m α̃p ∆α̃p Equidistant cylindrical 1 40.00 45.00 195 +0.0% 190 +0.0% 0 - 0 - 0 - 0 - 2 40.00 45.00 195 +0.0% 190 +0.0% 0 - 0 - 0 - 0 - 5 40.00 45.00 195 +0.0% 190 +0.0% 0 - 0 - 0 - 0 - Equidistant conic 1 45.00 0.68 195 +0.0% 10070 -0.2% 0 - 0.78 -21.2% 0 - 0.52 +0.0% 2 45.00 0.78 195 +0.0% 8778 +0.5% 0 - 0.90 -54.8% 0 - 0.60 +20.0% 5 45.00 2.74 195 +0.0% 2470 +0.0% 0 - 3.19 -27.7% 0 - 2.14 -0.9% Equidistant azimuthal 1 45.00 0.50 195 +0.0% 13718 -0.5% 0 - 0.50 -50.0% 0 - 0.50 +0.0% 2 45.00 0.71 195 +0.0% 9652 -0.6% 0 - 1.06 -79.3% 0 - 0.71 +0.0% 5 45.00 2.74 195 +0.0% 2508 +0.1% 0 - 4.11 -14.0% 0 - 2.75 -0.7% Bonne 1 0.42 0.99 16770 -4.1% 6916 -6.0% 8.53 +753.0% 1.14 +14.0% 0.38 +0.0% 0.44 +0.0% 2 0.85 1.59 8307 -0.7% 4332 +0.0% 9.96 +295.2% 1.83 -8.0% 0.76 +0.0% 0.71 +0.0% 5 2.05 4.23 3471 -0.3% 1672 +1.2% 9.61 +93.0% 3.24 -34.7% 1.79 -2.7% 1.87 +2.8% Hassler 1 0.36 0.70 19539 -1.8% 9842 +0.2% 0.82 -18.0% 0.70 -30.0% 0.35 +0.0% 0.46 +0.0% 2 0.85 1.04 8307 -0.9% 6650 -0.3% 1.93 -3.5% 1.04 -50.7% 0.83 +1.2% 0.68 +1.5% 5 1.73 3.43 4134 -1.4% 2014 +0.6% 3.92 -21.0% 5.14 +8.0% 1.69 +1.8% 2.25 -3.0% Nicolosi 1 0.45 1.17 15200 +0.2% 6200 -0.4% 1.13 +13.0% 0.87 -13.0% 0.36 -10.0% 0.38 -5.0% 2 0.85 2.51 8132 -0.3% 2920 +0.3% 1.42 -29.0% 1.87 -6.0% 0.75 +1.4% 0.98 +21.0% 5 2.00 5.40 3496 +0.2% 1360 +0.6% 3.34 -32.2% 5.96 +21.9% 1.78 +0.6% 1.75 -0.6% Littrow 1 0.77 0.45 8112 -0.1% 13718 -0.6% 4.94 +394.0% 2.88 +288.0% 0.40 +0.0% 0.47 +4.4% 2 1.66 0.86 4329 -1.0% 7980 +0.6% 9.39 +369.5% 4.95 +150.0% 0.76 +1.3% 0.86 +11.7% 5 3.28 2.38 2184 +1.5% 2926 +0.8% 9.49 +95.3% 9.54 +91.6% 1.54 -0.7% 2.38 +11.7% Adams I. 1 0.88 1.51 8073 -0.2% 4598 +0.8% 7.39 +739.0% 5.02 +412.2% 0.43 +13.2% 0.31 -22.5% 2 1.52 2.73 4719 +0.6% 2622 +0.1% 9.26 +363.0% 8.29 +316.6% 0.74 +10.5% 0.54 -25.0% 5 3.35 6.49 2145 -0.1% 1064 +0.6% 9.96 +103.7% 9.00 +102.7% 1.62 +8.7% 1.26 -26.7% Geoinformatics FCE CTU 17(2), 2018 53 T. Bayer: Plotting the map projection graticule with discontinuities a) b) Figure 5.1: Comparison of the uniform a) and combined b) sampling (S3) for the Bonne projection, α = 5◦; the adaptive polygonal approximation is smoother. the higher-curvature regions, the combined sampling technique provides a smoother approxi- mation 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 high- curvature regions (see Bonne projection in Fig. 5.1). In general, uniform sampling preserves the curvature worse and provides more redundant data. 5.2. 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. If (λ− λ0) = 0 or |ϕ| = π/2, X = 0, Y = Rϕ. If ϕ = 0, X = R (λ− λ0) , Y = 0. If |λ− λ0| = π/2, X = R (λ− λ0) cosϕ, Y = R π 2 sinϕ. Otherwise, C = π2 4 , P = π |sinϕ| , S = C − ϕ2 P − 2 |ϕ| , A = (λ− λ0)2 C − 1, B = S2 −A[C − PS − (λ− λ0)2], and Y = R signϕ √ B − S A , X = R(λ− λ0) √ 1− Y 2 R2C, where Y takes the sign of ϕ. The singularities are shown in Fig. 5.2. Geoinformatics FCE CTU 17(2), 2018 54 T. Bayer: Plotting the map projection graticule with discontinuities -1 80 .0 -1 70 .0 -1 60 .0 -1 50 .0 -1 40 .0 -1 30 .0 -1 20 .0 -1 10 .0 -1 00 .0 -9 0. 0 -8 0. 0 -7 0. 0 -6 0. 0 -5 0. 0 -4 0. 0 -3 0. 0 -2 0. 0 -1 0. 0 0. 0 10 .0 2 0 .0 30 .0 40 .0 5 0 .0 60 .0 70 .0 80 .0 90 .0 1 0 0. 0 1 1 0. 0 12 0. 0 13 0. 0 14 0. 0 15 0. 0 16 0. 0 17 0. 0 18 0. 0 -0 .0 0. 0 -90.0-90.0 -8 0.0 -80.0 -7 0.0 -70.0 -60 .0 -60.0 -50 .0 -50.0 -40 .0 -40.0 -30.0 -30.0 -20.0 -20.0 -10.0 -10.0 0.0 0.0 10.0 10.0 20.0 20.0 30.0 30.0 40.0 40.0 50.0 50. 0 60.0 60 .0 70.0 70 .0 80.0 80 .0 90.090.0 Figure 5.2: Fournier I. projection with the highlighted singularities leading to the subdivision into the tiles; normal aspect of the projection. Testing a normal aspect of the projection. Suppose that the graticule will be generated over the entire planisphere, so Ω = Ωϕ ×Ωλ, Ωϕ = [−π/2, π/2], Ωλ = [−π, π] with the offset of meridians and parallels ∆ϕ = ∆λ = 15◦, the sampling step is α = 1◦. Let us describe the selected steps of the algorithm; see Fig. 5.3: • Initialize “good data” as Ωg = Ω and push S ← Ωg. Pop the current interval Ωj = Ωϕ,j × Ωλ,j from S, detect the removable discontinuity at ϕ = 0, split Ωj to Ωϕ,j,1 = [−π/2, 0), Ωϕ,j,2 = (0, π/2]; see Step b). Create the new subintervalsΩj,1 = Ωϕ,j,1×Ωλ,j , Ω,j,2 = Ωϕ,j,2 ×Ωλ,j , and push S ← Ωj,1, S ← Ωj,2. • Due to the discontinuity at λ = −π/2, the current interval Ωj popped from the stack is split to Ωλ,j,1 = [−π,−π/2), Ωϕ,λ,2 = (−π/2, π]; see Step c). The created subintervals Ωj,1 = Ωϕ,j ×Ωλ,j,1, Ωj,2 = Ωϕ,j ×Ωλ,j,2 are pushed to the stack S ← Ωj,1, S ← Ωj,2. • Subsequently, the current interval Ωj = (0,−π/2] × [−π,−π/2) popped out from the stack is free from singularities, the graticule fragment over Ωj is constructed; see Step d). • Unfortunately, the next current interval Ωj = (0,−π/2]×(−π/2, π] contains a removable singularity at λ = π/2. It is split to Ωλ,j,1 = [−π/2, π/2), Ωϕ,λ,2 = (π/2, π], the created subintervals Ωj,1 = Ωϕ,j×Ωλ,j,1, Ω,j,2 = Ωϕ,j×Ωλ,j,2 are pushed to the stack S ← Ωj,1, S ← Ωj,2; see Step e). • Subsequently, the current interval Ωj = (0,−π/2]× [−π/2, π/2) popped from the stack is free from singularities, the graticule fragment over Ωj is constructed; see Step f). The lack of discontinuities refers to the next popped interval Ωj = (0,−π/2]× [π/2, π), over which the graticule is constructed; see Step g). Geoinformatics FCE CTU 17(2), 2018 55 T. Bayer: Plotting the map projection graticule with discontinuities a) b) c) d) e) f ) g) h) i) j) k) l) Figure 5.3: Reconstruction of the graticule from the tiles generated along the set of Ωg, Fournier I projection. Geoinformatics FCE CTU 17(2), 2018 56 T. Bayer: Plotting the map projection graticule with discontinuities 1 1 0 .0 1 2 0 .0 1 3 5 .0 1 5 0 .0 1 6 5 .0 1 8 0 .0 -7 0 .0 -6 0 .0 -4 5. 0 -3 0. 0 -15 .0 0.0 15.0 30.0 45 .0 60 .0 75 .0 90 .0 1 0 5 .0 1 1 0 .0 20.020.0 -1 8 0 .0 -1 6 5 .0 -1 5 0 .0 -1 3 5 .0 -1 2 0 .0 -1 0 5 .0 -9 0 .0 -7 5 .0 -7 0 .0 -1 6 0 .0 -1 6 0 .0 -90.0 -75.0 -60.0 -45.0 -30.0 -15.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 -90.0-90.0 -75.0 -75.0 -60 .0 -60.0 -4 5.0 -45.0 -3 0. 0 -30.0 -1 5 .0 -1 5 .0 0 .0 0 .0 1 5 .0 1 5 .0 30 .0 30.0 45 .0 45.0 60 .0 60.0 75 .0 75.0 90 .0 90.0 -90.0 -90. 0 -75.0 -75.0 -60.0 -60.0 -45.0 -45 .0 -30.0 -30 .0 -15.0 -15 .0 0.0 0.0 15.0 15 .0 30.0 30 .0 45.0 45 .0 60.0 60 .0 75.0 75 .0 90.090 .0 Figure 5.4: Fournier I. projection with the highlighted singularities leading to the subdivision into the tiles; an oblique aspect of the projection, K = [50◦, 20◦]. Testing an oblique aspect of the projection. For the oblique aspect of the projection, several new discontinuities appear. Suppose the arbitrary pole K = [50◦, 20◦], the central meridian shift is λ′0 = 0◦. Due to the latitude singularity at ϕk = 50◦, the latitude interval Ωϕ,j will be split into two parts: Ωϕ,j,1 = [−90◦, 50◦), Ωϕ,j,2 = (50◦, 90◦]. Analogously, the longitude singularity at λk − 180◦ = −160◦ splits the interval Ωλ,j into the three parts: Ωλ,j,1 = [−180◦,−160◦), Ωλ,j,2 = (−160◦, 15◦), Ωλ,j.3 = (15◦, 180◦]. The reconstructed graticule with the highlighted singularities is depicted in Fig. 5.4. 5.3. 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 ns, processed tiles nt, interval resizing nr, deleted intervals nd, and possible failures of the algorithm will be investigated. The impact of the central meridian shift λ′0. To verify the properties, behavior and numerical stability of the procedures, the following shifts λ′0 = {0◦, 30◦, ..., 150◦, 180◦} of the Geoinformatics FCE CTU 17(2), 2018 57 T. Bayer: Plotting the map projection graticule with discontinuities Table 3: The amount of splits ns, processed tiles nt, interval resizing nr, and deleted intervals nd for different values of λ′0 representing multiples of 30◦. Projection λ′0 0◦ 30◦ 60◦ 90◦ 120◦ 150◦ 180◦ Fournier I. nt 7 3 3 11 3 3 7 ns 3 1 1 5 1 1 3 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Van der Grinten I. nt 7 3 3 3 3 3 3 ns 3 1 1 1 1 1 1 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Ortelius Oval nt 3 1 1 1 1 1 1 ns 1 1 1 1 1 1 1 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Nicolosi nt 7 3 3 7 3 3 3 ns 3 1 1 3 1 1 1 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Hassler nt 1 1 1 1 1 1 1 ns 1 1 1 1 1 1 1 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Table 4: The amount of splits ns, processed tiles nt, interval resizing nr, and deleted intervals nd for different values of λ′0 which are not the exact multiples of 30◦. Projection λ′0 0.001◦ 30.001◦ 60.001◦ 90.001◦ 120.001◦ 150.001◦ 180.001◦ Fournier I. nt 5 3 3 8 3 3 5 ns 1 1 1 3 1 1 1 nr 1 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Van der Grinten I. nt 4 3 3 4 3 3 4 ns 1 1 1 1 1 1 1 nr 0 0 0 1 0 0 0 nd 0 0 0 0 0 0 0 Ortelius Oval nt 2 1 1 2 1 1 1 ns 0 0 0 0 0 0 0 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Nicolosi nt 5 3 3 8 3 3 5 ns 1 1 1 3 1 1 1 nr 1 0 0 0 0 0 1 nd 0 0 0 0 0 0 0 Hassler nt 4 3 3 3 3 3 4 ns 1 1 1 1 1 1 1 nr 0 0 0 0 0 0 0 nd 0 0 0 0 0 0 0 Geoinformatics FCE CTU 17(2), 2018 58 T. Bayer: Plotting the map projection graticule with discontinuities Table 5: The amount of splits ns, processed tiles nt, interval resizing nr, deleted intervals nd depending on λ′0 and the threshold ε. Projection λ′0 0.001◦ 30.001◦ 60.001◦ 90.001◦ 120.001◦ 150.001◦ 179.999◦ Fournier I. nt 4/17/145/ 1435 3/17/145/ 1435 3/17/145/ 1435 13/65/661/ 6591 3/17/145/ 1435 3/17/145/ 1435 4/17/145/ 1435 ns 1/1/1/2 1/1/1/2 1/1/1/2 5/5/5/6 1/1/1/2 1/1/1/2 1/1/1/2 nr 1/7/71/ 1212 0/7/71/ 1212 0/7/71/ 1212 1/27/325/ 5786 0/7/71/ 1212 0/7/71/ 1212 1/7/71/ 1212 nd 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 Van der Grinten I. nt 5/27/265/ 2633 5/27/259/ 2577 5/27/259/ 2577 5/27/259/ 2577 5/27/259/ 2577 5/27/259/ 2577 5/27/259/ 2633 ns 1/1/3/4 1/1/1/2 1/1/1/2 1/1/1/2 1/1/1/2 1/1/1/2 1/1/1/2 nr 1/12/130/ 2311 1/12/128/ 2283 1/12/128/ 2283 1/12/128/ 2283 1/12/128/ 2283 1/12/128/ 2283 1/12/128/ 2311 nd 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 Ortelius Oval nt 2/2/5/58 1/1/1/1 1/1/1/1 2/2/1/1 1/1/1/1 1/1/1/1 2/2/5/58 ns 0/0/1/3 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/1/3 nr 0/0/1/38 0/0/0/0 0/0/0/0 1/1/0/0 0/0/0/0 0/0/0/0 1/1/2/39 nd 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 Nicolosi nt 5/47/453/ 4519 3/47/453/ 4519 3/47/453/ 4519 13/95/969/ 9576 3/47/453/ 4519 3/47/453/ 4519 4/47/453/ 4519 ns 1/1/1/2 1/1/1/2 1/1/1/2 5/5/5/6 1/1/1/2 1/1/1/2 1/1/1/2 nr 1/22/225/ 4254 0/22/225/ 4254 0/22/225/ 4254 1/42/479/ 8828 0/22/225/ 4254 0/22/225/ 4254 1/22/225/ 4254 nd 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/0/1 0/0/1/1 Hassler nt 3/3/5/31 3/3/5/31 3/3/5/31 3/3/5/31 3/3/5/31 3/3/5/31 3/3/5/31 ns 1/1/1/1 1/1/1/1 1/1/1/1 1/1/1/1 1/1/1/1 1/1/1/1 1/1/1/1 nr 0/0/1/14 0/0/1/14 0/0/1/14 0/0/1/14 0/0/1/14 0/0/1/14 0/0/1/14 nd 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 0/0/0/0 central meridian representing multiples of 30◦, are involved. All projections are proposed in the oblique aspect, where K = [50◦, 20◦]. The results are summarized in Tab. 3. It is evident that the quantitative parameters of the algorithm depend on λ′0. If λ′0 coincides with the singularity, the number of interval corrections as well as the amount of processed tiles increase and vice versa; the maximum values can be found for λ′0 = {0◦, 90◦, 180◦}, which is in accordance with the detected singularities. Taking into account the complexity of equations, the less complicated formulas with few singularities are constructed from a smaller amount of tiles. Due to the position of K and λ′0, there are neither resized nor empty intervals. In general, only a few subdivisions are carried out. The impact of the numerical inaccuracy. The impact of numerical inaccuracies will also be illustrated. All projections are proposed in the transverse aspect, where K = [0◦, 90◦]. The central meridian shifts are not the exact multiples of 30◦, but λ′0 = {0.001◦, 30.001◦, ..., 150.001◦, 179.999◦}. Therefore, some resizing steps are expected. The results are summarized in Tab. 4. Comparing the quantitative parameters, the number of processed tiles nt decreased, but some intervals need to be resized. After this refinement they become correct; only a single resizing step is required. However, there was no need for a deletion of the empty interval. It is evident that “simple” equations are insensitive to the projection aspect. Geoinformatics FCE CTU 17(2), 2018 59 T. Bayer: Plotting the map projection graticule with discontinuities Figure 5.5: The list of reconstructed graticules of Fournier I., Van der Grinten I., Ortelius Oval, Nicolosi, and Hassler projections using combined sampling, α = 5◦; the central meridian shifts λ′0 = 90.001◦, λ′0 = 179.999◦. Geoinformatics FCE CTU 17(2), 2018 60 T. Bayer: Plotting the map projection graticule with discontinuities The impact of the central meridian shift λ′0 and the numerical threshold ε. The last test measures the quantitative parameters of the algorithms depending on the value of the threshold ε. This value has a strong influence on whether the result will be classified as a singularity. For testing, the following values are involved: ε = {1.0e−6, 1.0e−7, 1.0e−8, 1.0e−9}; the cartographic parameters remain unchained. The results are summarized in Tab. 5. It is evident that the amount of tiles nt strongly depends on ε, especially if a singularity at λ′ exists. Increasing ε by one order of magnitude, the amount of tiles nt changes in the same manner. Taking into account the results, it does not make a sense to use values ε < 1.0e−7. The algorithm is slowing down; its computational demands significantly grow up. It also becomes unstable and may lead to the stack overflow. Comparing the analyzed projections, the worst results are provided by the Nicolosi projection, where nt almost doesn’t depend on λ′0. In the most pessimistic case, almost 10 000 tiles are required. On the contrary, the Hassler projection with nt < 35 needs less effort. Summarize the facts, the reconstruction algorithm depends on the number of singularities as well as on the threshold ε. For ε ≥ 1.0e−6, only a few splits are required. All the reconstructed graticules of the analyzed projections for λ′0 = 90.001◦, and λ′0 = 179.999◦ can be found in Fig. 5.5. 6. Conclusion This article presented a new algorithm for the combined sampling of the projection graticule handling all constant values of the projection. The proposed method combines the uniform and adaptive sampling techniques with the triple recursive subdivision. It controls handling the discontinuities and their directions; the detection criterion is presented in Sec. 3.1. In general, the splitting procedure is sensitive to ε(if λ′0 → ε) as well as to the numerical errors of intersections m′(λ) and m(λ), m′(λ) and p(ϕ). From the topological point of view, the meridian and parallel intersections lack the vertex, which can cause some problems in GIS processing; the dissolve operation needs to be applied before. The discontinuity detection without its classification (or, a wrong classification) leads to the duplication of lines along the singularities. The idea of splitting the interval to the set of disjoint intervals without the singularities is efficient and works for any common projection and its coordinate functions. Compared to the uniform sampling, our experiments confirmed that the combined sampling techniques S1-S3 require fewer data in regions with the high curvature (more than five times lower values of α). For the circular arcs, circles or ellipses the uniform sampling is slightly more efficient. The lowest data redundancy refers to the straight lines, only 1% of the uniform sampling data is sufficient for the polygonal approximation. This is beneficial for the cylindrical, conic, azimuthal and pseudocylindrical projections. A minor disadvantage of the proposed solution is a slightly complicated implementation of the algorithms, the source code in C++/Java is available in the GitHub repository https://github.com/bayertom/graticule. Geoinformatics FCE CTU 17(2), 2018 61 https://github.com/bayertom/graticule T. Bayer: Plotting the map projection graticule with discontinuities References [1] Giampietro Allasia, Renata Besenghi, and Roberto Cavoretto. “Adaptive detection and approximation of unknown surface discontinuities from scattered data”. In: Simulation Modelling Practice and Theory 17.6 (2009), pp. 1059–1070. [2] Francesc Arandiga et al. “Interpolation and Approximation of Piecewise Smooth Func- tions”. In: SIAM Journal on Numerical Analysis 43.1 (2005), pp. 41–57. doi: 10.1137/ S0036142903426245. [3] Rick Archibald, Anne Gelb, and Jungho Yoon. “Determining the locations and disconti- nuities in the derivatives of functions”. In: Applied Numerical Mathematics 58.5 (2008), pp. 577–592. [4] Tomas Bayer. “Advanced methods for the estimation of an unknown projection from a map”. In:GeoInformatica 20.2 (2016), pp. 241–284. doi: 10.1007/s10707-015-0234-x. [5] Tomas Bayer. “Efficient plotting of functions with discontinuities based on combined sampling”. English. In: FCE GeoInformatics 18.3 (2018), pp. 621–669. issn: 1384-6175. doi: 10.1007/s10707-013-0200-4. [6] Tomas Bayer. “Estimation of an unknown cartographic projection and its parameters from the map”. English. In: GeoInformatica 18.3 (2014), pp. 621–669. issn: 1384-6175. doi: 10.1007/s10707-013-0200-4. [7] P. Binev et al. “Adaptive approximation of curves”. In: Approximation Theory (2004), pp. 43–57. [8] Mira Bozzini, Licia Lenarduzzi, and Milvia Rossini. “Non-regular Surface Approxima- tion”. In: Mathematical Methods for Curves and Surfaces: 8th International Conference, MMCS 2012, Oslo, Norway, June 28 – July 3, 2012, Revised Selected Papers. Ed. by Michael Floater et al. Berlin, Heidelberg: Springer Berlin Heidelberg, 2014, pp. 68–87. isbn: 978-3-642-54382-1. doi: 10.1007/978-3-642-54382-1_5. [9] Mira Bozzini, Licia Lenarduzzi, and Robert Schaback. “Adaptive Interpolation by Scaled Multiquadrics”. In: Advances in Computational Mathematics 16.4 (May 2002), pp. 375– 387. issn: 1572-9044. doi: 10.1023/A:1014584220418. [10] Filipe de Carvalho Nascimento et al. “Approximating implicit curves on plane and surface triangulations with affine arithmetic”. In: Computers & Graphics 40 (2014), pp. 36–48. [11] Andrew Crampton and John C Mason. “Detecting and approximating fault lines from randomly scattered data”. In: Numerical Algorithms 39.1 (2005), pp. 115–130. [12] Luiz Henrique de Figueiredo. “Adaptive sampling of parametric curves”. In: Graphics Gems V 5 (1995), pp. 173–178. [13] Tim Gutzmer and Armin Iske. “Detection of discontinuities in scattered data approxi- mation”. In: Numerical Algorithms 16.2 (1997), pp. 155–170. [14] John Krumm. Intersection of Two Planes. 2016. [15] David Lee. “Detection, Classification, and Measurement of Discontinuities”. In: SIAM Journal on Scientific and Statistical Computing 12.2 (1991), pp. 311–341. doi: 10. 1137/0912018. Geoinformatics FCE CTU 17(2), 2018 62 https://doi.org/10.1137/S0036142903426245 https://doi.org/10.1137/S0036142903426245 https://doi.org/10.1007/s10707-015-0234-x https://doi.org/10.1007/s10707-013-0200-4 https://doi.org/10.1007/s10707-013-0200-4 https://doi.org/10.1007/978-3-642-54382-1_5 https://doi.org/10.1023/A:1014584220418 https://doi.org/10.1137/0912018 https://doi.org/10.1137/0912018 T. Bayer: Plotting the map projection graticule with discontinuities [16] Licia Lenarduzzi and Robert Schaback. “Kernel-based adaptive approximation of func- tions with discontinuities”. In: Applied Mathematics and Computation 307.Supplement C (2017), pp. 113–123. issn: 0096-3003. doi: https://doi.org/10.1016/j.amc.2017. 02.043. [17] Hélio Lopes, João Batista Oliveira, and Luiz Henrique de Figueiredo. “Robust adaptive polygonal approximation of implicit curves”. In: Computers & Graphics 26.6 (2002), pp. 841–852. [18] María Cruz López de Silanes, María Cruz Parra, and Juan José Torrens. “On a New Characterization of Finite Jump Discontinuities and Its Application to Vertical Fault Detection”. In: Math. Comput. Simul. 77.2-3 (Mar. 2008), pp. 247–256. issn: 0378-4754. doi: 10.1016/j.matcom.2007.08.008. [19] María Cruz López de Silanes, María Cruz Parra, and Juan José Torrens. “Vertical and oblique fault detection in explicit surfaces”. In: Journal of Computational and Applied Mathematics 140.1 (2002). Int. Congress on Computational and Applied Mathematics 2000, pp. 559–585. issn: 0377-0427. doi: https://doi.org/10.1016/S0377-0427(01) 00601-X. [20] M Oliveria et al. “Universal high order subroutine with new shock detector for shock boundary layer interaction”. In: In other words 10 (2009), pp. 1–2. [21] Afonso Paiva et al. “Approximating implicit curves on triangulations with affine arith- metic”. In: Graphics, Patterns and Images (SIBGRAPI), 2012 25th SIBGRAPI Con- ference on. IEEE. 2012, pp. 94–101. [22] M. Rossini. “2D-discontinuity detection from scattered data”. In: Computing 61.3 (Sept. 1998), pp. 215–234. issn: 1436-5057. doi: 10.1007/BF02684351. [23] Moshe Shpitalni, Yoram Koren, and CC Lo. “Realtime curve interpolators”. In: Computer- Aided Design 26.11 (1994), pp. 832–838. [24] John P. Snyder. Map projections – A working manual. Tech. rep. 1395. U. S. Geological Survey, 1987. url: http://pubs.er.usgs.gov/usgspubs/pp/pp1395. [25] Daniel C.H. Yang and Tom Kong. “Parametric interpolator versus linear interpolator for precision CNC machining”. In: Computer-Aided Design 26.3 (1994). Special Issue:NC machining and cutter-path generation, pp. 225–234. issn: 0010-4485. doi: https:// doi.org/10.1016/0010-4485(94)90045-0. [26] S-S Yeh and P-L Hsu. “The speed-controlled interpolator for machining parametric curves”. In: Computer-Aided Design 31.5 (1999), pp. 349–357. [27] Syh-Shiuh Yeh and Pau-Lo Hsu. “Adaptive-feedrate interpolation for parametric curves with a confined chord error”. In: Computer-aided design 34.3 (2002), pp. 229–237. Geoinformatics FCE CTU 17(2), 2018 63 https://doi.org/https://doi.org/10.1016/j.amc.2017.02.043 https://doi.org/https://doi.org/10.1016/j.amc.2017.02.043 https://doi.org/10.1016/j.matcom.2007.08.008 https://doi.org/https://doi.org/10.1016/S0377-0427(01)00601-X https://doi.org/https://doi.org/10.1016/S0377-0427(01)00601-X https://doi.org/10.1007/BF02684351 http://pubs.er.usgs.gov/usgspubs/pp/pp1395 https://doi.org/https://doi.org/10.1016/0010-4485(94)90045-0 https://doi.org/https://doi.org/10.1016/0010-4485(94)90045-0 Geoinformatics FCE CTU 17(2), 2018 64 T. Bayer: Plotting the map projection graticule with discontinuities Introduction Related Work Map projection and its properties Map projection singularities Intersection of the great circle and meridian/parallel arcs Graticule reconstruction using the combined sampling Graticule construction with singularities Combined sampling of meridian and parallel fragments Experiments and results Comparison of uniform and combined sampling Detection of discontinuities: Fournier I. projection Measuring the quantitative parameters of the algorithm Conclusion