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