1"""
2    kartograph - a svg mapping library
3    Copyright (C) 2011  Gregor Aisch
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as
7    published by the Free Software Foundation, either version 3 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <http://www.gnu.org/licenses/>.
17"""
18
19from cylindrical import Cylindrical
20import math
21from math import radians as rad
22
23
24class PseudoCylindrical(Cylindrical):
25    def __init__(self, lon0=0.0, flip=0):
26        Cylindrical.__init__(self, lon0=lon0, flip=flip)
27
28
29class NaturalEarth(PseudoCylindrical):
30    """
31    src: http://www.shadedrelief.com/NE_proj/
32    """
33    def __init__(self, lat0=0.0, lon0=0.0, flip=0):
34        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
35        from math import pi
36        s = self
37        s.A0 = 0.8707
38        s.A1 = -0.131979
39        s.A2 = -0.013791
40        s.A3 = 0.003971
41        s.A4 = -0.001529
42        s.B0 = 1.007226
43        s.B1 = 0.015085
44        s.B2 = -0.044475
45        s.B3 = 0.028874
46        s.B4 = -0.005916
47        s.C0 = s.B0
48        s.C1 = 3 * s.B1
49        s.C2 = 7 * s.B2
50        s.C3 = 9 * s.B3
51        s.C4 = 11 * s.B4
52        s.EPS = 1e-11
53        s.MAX_Y = 0.8707 * 0.52 * pi
54
55    def project(self, lon, lat):
56        lon, lat = self.ll(lon, lat)
57        lplam = rad(lon)
58        lpphi = rad(lat * -1)
59        phi2 = lpphi * lpphi
60        phi4 = phi2 * phi2
61        x = lplam * (self.A0 + phi2 * (self.A1 + phi2 * (self.A2 + phi4 * phi2 * (self.A3 + phi2 * self.A4)))) * 180 + 500
62        y = lpphi * (self.B0 + phi2 * (self.B1 + phi4 * (self.B2 + self.B3 * phi2 + self.B4 * phi4))) * 180 + 270
63        return (x, y)
64
65
66class Robinson(PseudoCylindrical):
67
68    def __init__(self, lat0=0.0, lon0=0.0, flip=0):
69        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
70        self.X = [1, -5.67239e-12, -7.15511e-05, 3.11028e-06,  0.9986, -0.000482241, -2.4897e-05, -1.33094e-06, 0.9954, -0.000831031, -4.4861e-05, -9.86588e-07, 0.99, -0.00135363, -5.96598e-05, 3.67749e-06, 0.9822, -0.00167442, -4.4975e-06, -5.72394e-06, 0.973, -0.00214869, -9.03565e-05, 1.88767e-08, 0.96, -0.00305084, -9.00732e-05, 1.64869e-06, 0.9427, -0.00382792, -6.53428e-05, -2.61493e-06, 0.9216, -0.00467747, -0.000104566, 4.8122e-06, 0.8962, -0.00536222, -3.23834e-05, -5.43445e-06, 0.8679, -0.00609364, -0.0001139, 3.32521e-06, 0.835, -0.00698325, -6.40219e-05, 9.34582e-07, 0.7986, -0.00755337, -5.00038e-05, 9.35532e-07, 0.7597, -0.00798325, -3.59716e-05, -2.27604e-06, 0.7186, -0.00851366, -7.0112e-05, -8.63072e-06, 0.6732, -0.00986209, -0.000199572, 1.91978e-05, 0.6213, -0.010418, 8.83948e-05, 6.24031e-06, 0.5722, -0.00906601, 0.000181999, 6.24033e-06, 0.5322,  0.,  0.,  0.]
71        self.Y = [0, 0.0124, 3.72529e-10, 1.15484e-09, 0.062, 0.0124001, 1.76951e-08, -5.92321e-09, 0.124, 0.0123998, -7.09668e-08, 2.25753e-08, 0.186, 0.0124008, 2.66917e-07, -8.44523e-08, 0.248, 0.0123971, -9.99682e-07, 3.15569e-07, 0.31, 0.0124108, 3.73349e-06, -1.1779e-06, 0.372, 0.0123598, -1.3935e-05, 4.39588e-06, 0.434, 0.0125501, 5.20034e-05, -1.00051e-05, 0.4968, 0.0123198, -9.80735e-05, 9.22397e-06, 0.5571, 0.0120308, 4.02857e-05, -5.2901e-06, 0.6176, 0.0120369, -3.90662e-05, 7.36117e-07, 0.6769, 0.0117015, -2.80246e-05, -8.54283e-07, 0.7346, 0.0113572, -4.08389e-05, -5.18524e-07, 0.7903, 0.0109099, -4.86169e-05, -1.0718e-06, 0.8435, 0.0103433, -6.46934e-05, 5.36384e-09, 0.8936, 0.00969679, -6.46129e-05, -8.54894e-06, 0.9394, 0.00840949, -0.000192847, -4.21023e-06, 0.9761, 0.00616525, -0.000256001, -4.21021e-06, 1.,  0.,  0.,  0]
72        self.NODES = 18
73        self.FXC = 0.8487
74        self.FYC = 1.3523
75        self.C1 = 11.45915590261646417544
76        self.RC1 = 0.08726646259971647884
77        self.ONEEPS = 1.000001
78        self.EPS = 1e-8
79
80    def _poly(self, arr, off, z):
81        return arr[off] + z * (arr[off + 1] + z * (arr[off + 2] + z * (arr[off + 3])))
82
83    def project(self, lon, lat):
84        lon, lat = self.ll(lon, lat)
85        lplam = rad(lon)
86        lpphi = rad(lat * -1)
87
88        phi = abs(lpphi)
89        i = int(phi * self.C1)
90        if i >= self.NODES:
91            i = self.NODES - 1
92        phi = math.degrees(phi - self.RC1 * i)
93        i *= 4
94        x = 1000 * self._poly(self.X, i, phi) * self.FXC * lplam
95        y = 1000 * self._poly(self.Y, i, phi) * self.FYC
96        if lpphi < 0.0:
97            y = -y
98
99        return (x, y)
100
101
102class EckertIV(PseudoCylindrical):
103
104    def __init__(self, lon0=0.0, lat0=0, flip=0):
105        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
106
107        self.C_x = .42223820031577120149
108        self.C_y = 1.32650042817700232218
109        self.RC_y = .75386330736002178205
110        self.C_p = 3.57079632679489661922
111        self.RC_p = .28004957675577868795
112        self.EPS = 1e-7
113        self.NITER = 6
114
115    def project(self, lon, lat):
116        lon, lat = self.ll(lon, lat)
117        lplam = rad(lon)
118        lpphi = rad(lat * -1)
119
120        p = self.C_p * math.sin(lpphi)
121        V = lpphi * lpphi
122        lpphi *= 0.895168 + V * (0.0218849 + V * 0.00826809)
123
124        i = self.NITER
125        while i > 0:
126            c = math.cos(lpphi)
127            s = math.sin(lpphi)
128            V = (lpphi + s * (c + 2.) - p) / (1. + c * (c + 2.) - s * s)
129            lpphi -= V
130            if abs(V) < self.EPS:
131                break
132            i -= 1
133
134        if i == 0:
135            x = self.C_x * lplam
136            y = (self.C_y, - self.C_y)[lpphi < 0]
137        else:
138            x = self.C_x * lplam * (1. + math.cos(lpphi))
139            y = self.C_y * math.sin(lpphi)
140        return (x, y)
141
142
143class Sinusoidal(PseudoCylindrical):
144
145    def __init__(self, lon0=0.0, lat0=0.0, flip=0):
146        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
147
148    def project(self, lon, lat):
149        lon, lat = self.ll(lon, lat)
150        lam = rad(lon)
151        phi = rad(lat * -1)
152        x = 1032 * lam * math.cos(phi)
153        y = 1032 * phi
154        return (x, y)
155
156
157class Mollweide(PseudoCylindrical):
158
159    def __init__(self, p=1.5707963267948966, lon0=0.0, lat0=0.0, cx=None, cy=None, cp=None, flip=0):
160        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
161        self.MAX_ITER = 10
162        self.TOLERANCE = 1e-7
163
164        if p != None:
165            p2 = p + p
166            sp = math.sin(p)
167            r = math.sqrt(math.pi * 2.0 * sp / (p2 + math.sin(p2)))
168            self.cx = 2. * r / math.pi
169            self.cy = r / sp
170            self.cp = p2 + math.sin(p2)
171        elif cx != None and cy != None and cp != None:
172            self.cx = cx
173            self.cy = cy
174            self.cp = cp
175        else:
176            assert False, 'either p or cx,cy,cp must be defined'
177
178    def project(self, lon, lat):
179        lon, lat = self.ll(lon, lat)
180        lam = rad(lon)
181        phi = rad(lat)
182
183        k = self.cp * math.sin(phi)
184        i = self.MAX_ITER
185
186        while i != 0:
187            v = (phi + math.sin(phi) - k) / (1. + math.cos(phi))
188            phi -= v
189            if abs(v) < self.TOLERANCE:
190                break
191            i -= 1
192
193        if i == 0:
194            phi = (self.HALFPI, -self.HALFPI)[phi < 0]
195        else:
196            phi *= 0.5
197
198        x = 1000 * self.cx * lam * math.cos(phi)
199        y = 1000 * self.cy * math.sin(phi)
200        return (x, y * -1)
201
202
203class GoodeHomolosine(PseudoCylindrical):
204
205    def __init__(self, lon0=0, flip=0):
206        self.lat1 = 41.737
207        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
208        self.p1 = Mollweide()
209        self.p0 = Sinusoidal()
210
211    def project(self, lon, lat):
212        lon, lat = self.ll(lon, lat)
213        #lon = me.clon(lon)
214        if abs(lat) > self.lat1:
215            return self.p1.project(lon, lat)
216        else:
217            return self.p0.project(lon, lat)
218
219
220class WagnerIV(Mollweide):
221    def __init__(self, lon0=0, lat0=0, flip=0):
222        # p=math.pi/3
223        Mollweide.__init__(self, p=1.0471975511965976, flip=flip)
224
225
226class WagnerV(Mollweide):
227    def __init__(self, lat0=0, lon0=0, flip=0):
228        Mollweide.__init__(self, cx=0.90977, cy=1.65014, cp=3.00896, flip=flip)
229
230
231class Loximuthal(PseudoCylindrical):
232
233    minLat = -89
234    maxLat = 89
235
236    def __init__(self, lon0=0.0, lat0=0.0, flip=0):
237        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
238        if flip == 1:
239            lat0 = -lat0
240        self.lat0 = lat0
241        self.phi0 = rad(lat0)
242
243    def project(self, lon, lat):
244        lon, lat = self.ll(lon, lat)
245        lam = rad(lon)
246        phi = rad(lat)
247        if phi == self.phi0:
248            x = lam * math.cos(self.phi0)
249        else:
250            try:
251                x = lam * (phi - self.phi0) / (math.log(math.tan(self.QUARTERPI + phi * 0.5)) - math.log(math.tan(self.QUARTERPI + self.phi0 * 0.5)))
252            except:
253                return None
254        x *= 1000
255        y = 1000 * (phi - self.phi0)
256        return (x, y * -1)
257
258    def attrs(self):
259        p = super(Loximuthal, self).attrs()
260        p['lat0'] = self.lat0
261        return p
262
263    @staticmethod
264    def attributes():
265        return ['lon0', 'lat0', 'flip']
266
267
268class CantersModifiedSinusoidalI(PseudoCylindrical):
269    """
270    Canters, F. (2002) Small-scale Map projection Design. p. 218-219.
271    Modified Sinusoidal, equal-area.
272
273    implementation borrowed from
274    http://cartography.oregonstate.edu/temp/AdaptiveProjection/src/projections/Canters1.js
275    """
276
277    def __init__(self, lon0=0.0, flip=0):
278        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
279        self.C1 = 1.1966
280        self.C3 = -0.1290
281        self.C3x3 = 3 * self.C3
282        self.C5 = -0.0076
283        self.C5x5 = 5 * self.C5
284
285    def project(self, lon, lat):
286        me = self
287        lon, lat = me.ll(lon, lat)
288
289        lon = rad(lon)
290        lat = rad(lat)
291
292        y2 = lat * lat
293        y4 = y2 * y2
294        x = 1000 * lon * math.cos(lat) / (me.C1 + me.C3x3 * y2 + me.C5x5 * y4)
295        y = 1000 * lat * (me.C1 + me.C3 * y2 + me.C5 * y4)
296        return (x, y * -1)
297
298
299class Hatano(PseudoCylindrical):
300
301    def __init__(me, lon0=0, flip=0):
302        PseudoCylindrical.__init__(me, lon0=lon0, flip=flip)
303        me.NITER = 20
304        me.EPS = 1e-7
305        me.ONETOL = 1.000001
306        me.CN = 2.67595
307        me.CS = 2.43763
308        me.RCN = 0.37369906014686373063
309        me.RCS = 0.41023453108141924738
310        me.FYCN = 1.75859
311        me.FYCS = 1.93052
312        me.RYCN = 0.56863737426006061674
313        me.RYCS = 0.51799515156538134803
314        me.FXC = 0.85
315        me.RXC = 1.17647058823529411764
316
317    def project(me, lon, lat):
318        [lon, lat] = me.ll(lon, lat)
319        lam = rad(lon)
320        phi = rad(lat)
321        c = math.sin(phi) * (me.CN, me.CS)[phi < 0.0]
322        for i in range(me.NITER, 0, -1):
323            th1 = (phi + math.sin(phi) - c) / (1.0 + math.cos(phi))
324            phi -= th1
325            if abs(th1) < me.EPS:
326                break
327        phi *= 0.5
328        x = 1000 * me.FXC * lam * math.cos(phi)
329        y = 1000 * math.sin(phi) * (me.FYCN, me.FYCS)[phi < 0.0]
330        return (x, y * -1)
331
332
333class Aitoff(PseudoCylindrical):
334    """
335    Aitoff projection
336
337    implementation taken from
338    Snyder, Map projections - A working manual
339    """
340    def __init__(self, lon0=0, flip=0):
341        PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
342        self.winkel = False
343        self.COSPHI1 = 0.636619772367581343
344
345    def project(me, lon, lat):
346        [lon, lat] = me.ll(lon, lat)
347        lam = rad(lon)
348        phi = rad(lat)
349        c = 0.5 * lam
350        d = math.acos(math.cos(phi) * math.cos(c))
351        if d != 0:
352            y = 1.0 / math.sin(d)
353            x = 2.0 * d * math.cos(phi) * math.sin(c) * y
354            y *= d * math.sin(phi)
355        else:
356            x = y = 0
357        if me.winkel:
358            x = (x + lam * me.COSPHI1) * 0.5
359            y = (y + phi) * 0.5
360        return (x * 1000, y * -1000)
361
362
363class Winkel3(Aitoff):
364
365    def __init__(self, lon0=0, flip=0):
366        Aitoff.__init__(self, lon0=lon0, flip=flip)
367        self.winkel = True
368
369
370class Nicolosi(PseudoCylindrical):
371
372    def __init__(me, lon0=0, flip=0):
373        me.EPS = 1e-10
374        PseudoCylindrical.__init__(me, lon0=lon0, flip=flip)
375        me.r = me.HALFPI * 100
376        sea = []
377        r = me.r
378        for phi in range(0, 361):
379            sea.append((math.cos(rad(phi)) * r, math.sin(rad(phi)) * r))
380        me.sea = sea
381
382    def _clon(me, lon):
383        lon -= me.lon0
384        if lon < -180:
385            lon += 360
386        elif lon > 180:
387            lon -= 360
388        return lon
389
390    def _visible(me, lon, lat):
391        #lon = me._clon(lon)
392        return lon > -90 and lon < 90
393
394    def _truncate(me, x, y):
395        theta = math.atan2(y, x)
396        x1 = me.r * math.cos(theta)
397        y1 = me.r * math.sin(theta)
398        return (x1, y1)
399
400    def world_bounds(self, bbox, llbbox=(-180, -90, 180, 90)):
401        if llbbox == (-180, -90, 180, 90):
402            d = self.r * 2
403            bbox.update((-d, -d))
404            bbox.update((d, d))
405        else:
406            bbox = super(PseudoCylindrical, self).world_bounds(bbox, llbbox)
407        return bbox
408
409    def sea_shape(self, llbbox=(-180, -90, 180, 90)):
410        out = []
411        if llbbox == (-180, -90, 180, 90) or llbbox == [-180, -90, 180, 90]:
412            for phi in range(0, 360):
413                x = math.cos(math.radians(phi)) * self.r
414                y = math.sin(math.radians(phi)) * self.r
415                out.append((x, y))
416            out = [out]
417        else:
418            out = super(PseudoCylindrical, self).sea_shape(llbbox)
419        return out
420
421    def project(me, lon, lat):
422        [lon, lat] = me.ll(lon, lat)
423        lam = rad(lon)
424        phi = rad(lat)
425
426        if abs(lam) < me.EPS:
427            x = 0
428            y = phi
429        elif abs(phi) < me.EPS:
430            x = lam
431            y = 0
432        elif abs(abs(lam) - me.HALFPI) < me.EPS:
433            x = lam * math.cos(phi)
434            y = me.HALFPI * math.sin(phi)
435        elif abs(abs(phi) - me.HALFPI) < me.EPS:
436            x = 0
437            y = phi
438        else:
439            tb = me.HALFPI / lam - lam / me.HALFPI
440            c = phi / me.HALFPI
441            sp = math.sin(phi)
442            d = (1 - c * c) / (sp - c)
443            r2 = tb / d
444            r2 *= r2
445            m = (tb * sp / d - 0.5 * tb) / (1.0 + r2)
446            n = (sp / r2 + 0.5 * d) / (1.0 + 1.0 / r2)
447            x = math.cos(phi)
448            x = math.sqrt(m * m + x * x / (1.0 + r2))
449            x = me.HALFPI * (m + (x, -x)[lam < 0])
450            f = n * n - (sp * sp / r2 + d * sp - 1.0) / (1.0 + 1.0 / r2)
451            if f < 0:
452                y = phi
453            else:
454                y = math.sqrt(f)
455                y = me.HALFPI * (n + (-y, y)[phi < 0])
456        return (x * 100, y * -100)
457
458    def plot(self, polygon, truncate=True):
459        polygons = self._shift_polygon(polygon)
460        plotted = []
461        for polygon in polygons:
462            points = []
463            ignore = True
464            for (lon, lat) in polygon:
465                vis = self._visible(lon, lat)
466                if vis:
467                    ignore = False
468                x, y = self.project(lon, lat)
469                if not vis and truncate:
470                    points.append(self._truncate(x, y))
471                else:
472                    points.append((x, y))
473            if ignore:
474                continue
475            plotted.append(points)
476        return plotted
477