Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np 

2 

3from matplotlib import cbook, rcParams 

4from matplotlib.axes import Axes 

5import matplotlib.axis as maxis 

6from matplotlib.patches import Circle 

7from matplotlib.path import Path 

8import matplotlib.spines as mspines 

9from matplotlib.ticker import ( 

10 Formatter, NullLocator, FixedLocator, NullFormatter) 

11from matplotlib.transforms import Affine2D, BboxTransformTo, Transform 

12 

13 

14class GeoAxes(Axes): 

15 """An abstract base class for geographic projections.""" 

16 class ThetaFormatter(Formatter): 

17 """ 

18 Used to format the theta tick labels. Converts the native 

19 unit of radians into degrees and adds a degree symbol. 

20 """ 

21 def __init__(self, round_to=1.0): 

22 self._round_to = round_to 

23 

24 def __call__(self, x, pos=None): 

25 degrees = round(np.rad2deg(x) / self._round_to) * self._round_to 

26 if rcParams['text.usetex'] and not rcParams['text.latex.unicode']: 

27 return r"$%0.0f^\circ$" % degrees 

28 else: 

29 return "%0.0f\N{DEGREE SIGN}" % degrees 

30 

31 RESOLUTION = 75 

32 

33 def _init_axis(self): 

34 self.xaxis = maxis.XAxis(self) 

35 self.yaxis = maxis.YAxis(self) 

36 # Do not register xaxis or yaxis with spines -- as done in 

37 # Axes._init_axis() -- until GeoAxes.xaxis.cla() works. 

38 # self.spines['geo'].register_axis(self.yaxis) 

39 self._update_transScale() 

40 

41 def cla(self): 

42 Axes.cla(self) 

43 

44 self.set_longitude_grid(30) 

45 self.set_latitude_grid(15) 

46 self.set_longitude_grid_ends(75) 

47 self.xaxis.set_minor_locator(NullLocator()) 

48 self.yaxis.set_minor_locator(NullLocator()) 

49 self.xaxis.set_ticks_position('none') 

50 self.yaxis.set_ticks_position('none') 

51 self.yaxis.set_tick_params(label1On=True) 

52 # Why do we need to turn on yaxis tick labels, but 

53 # xaxis tick labels are already on? 

54 

55 self.grid(rcParams['axes.grid']) 

56 

57 Axes.set_xlim(self, -np.pi, np.pi) 

58 Axes.set_ylim(self, -np.pi / 2.0, np.pi / 2.0) 

59 

60 def _set_lim_and_transforms(self): 

61 # A (possibly non-linear) projection on the (already scaled) data 

62 self.transProjection = self._get_core_transform(self.RESOLUTION) 

63 

64 self.transAffine = self._get_affine_transform() 

65 

66 self.transAxes = BboxTransformTo(self.bbox) 

67 

68 # The complete data transformation stack -- from data all the 

69 # way to display coordinates 

70 self.transData = \ 

71 self.transProjection + \ 

72 self.transAffine + \ 

73 self.transAxes 

74 

75 # This is the transform for longitude ticks. 

76 self._xaxis_pretransform = \ 

77 Affine2D() \ 

78 .scale(1, self._longitude_cap * 2) \ 

79 .translate(0, -self._longitude_cap) 

80 self._xaxis_transform = \ 

81 self._xaxis_pretransform + \ 

82 self.transData 

83 self._xaxis_text1_transform = \ 

84 Affine2D().scale(1, 0) + \ 

85 self.transData + \ 

86 Affine2D().translate(0, 4) 

87 self._xaxis_text2_transform = \ 

88 Affine2D().scale(1, 0) + \ 

89 self.transData + \ 

90 Affine2D().translate(0, -4) 

91 

92 # This is the transform for latitude ticks. 

93 yaxis_stretch = Affine2D().scale(np.pi * 2, 1).translate(-np.pi, 0) 

94 yaxis_space = Affine2D().scale(1, 1.1) 

95 self._yaxis_transform = \ 

96 yaxis_stretch + \ 

97 self.transData 

98 yaxis_text_base = \ 

99 yaxis_stretch + \ 

100 self.transProjection + \ 

101 (yaxis_space + 

102 self.transAffine + 

103 self.transAxes) 

104 self._yaxis_text1_transform = \ 

105 yaxis_text_base + \ 

106 Affine2D().translate(-8, 0) 

107 self._yaxis_text2_transform = \ 

108 yaxis_text_base + \ 

109 Affine2D().translate(8, 0) 

110 

111 def _get_affine_transform(self): 

112 transform = self._get_core_transform(1) 

113 xscale, _ = transform.transform((np.pi, 0)) 

114 _, yscale = transform.transform((0, np.pi/2)) 

115 return Affine2D() \ 

116 .scale(0.5 / xscale, 0.5 / yscale) \ 

117 .translate(0.5, 0.5) 

118 

119 def get_xaxis_transform(self, which='grid'): 

120 cbook._check_in_list(['tick1', 'tick2', 'grid'], which=which) 

121 return self._xaxis_transform 

122 

123 def get_xaxis_text1_transform(self, pad): 

124 return self._xaxis_text1_transform, 'bottom', 'center' 

125 

126 def get_xaxis_text2_transform(self, pad): 

127 return self._xaxis_text2_transform, 'top', 'center' 

128 

129 def get_yaxis_transform(self, which='grid'): 

130 cbook._check_in_list(['tick1', 'tick2', 'grid'], which=which) 

131 return self._yaxis_transform 

132 

133 def get_yaxis_text1_transform(self, pad): 

134 return self._yaxis_text1_transform, 'center', 'right' 

135 

136 def get_yaxis_text2_transform(self, pad): 

137 return self._yaxis_text2_transform, 'center', 'left' 

138 

139 def _gen_axes_patch(self): 

140 return Circle((0.5, 0.5), 0.5) 

141 

142 def _gen_axes_spines(self): 

143 return {'geo': mspines.Spine.circular_spine(self, (0.5, 0.5), 0.5)} 

144 

145 def set_yscale(self, *args, **kwargs): 

146 if args[0] != 'linear': 

147 raise NotImplementedError 

148 

149 set_xscale = set_yscale 

150 

151 def set_xlim(self, *args, **kwargs): 

152 raise TypeError("It is not possible to change axes limits " 

153 "for geographic projections. Please consider " 

154 "using Basemap or Cartopy.") 

155 

156 set_ylim = set_xlim 

157 

158 def format_coord(self, lon, lat): 

159 'return a format string formatting the coordinate' 

160 lon, lat = np.rad2deg([lon, lat]) 

161 if lat >= 0.0: 

162 ns = 'N' 

163 else: 

164 ns = 'S' 

165 if lon >= 0.0: 

166 ew = 'E' 

167 else: 

168 ew = 'W' 

169 return ('%f\N{DEGREE SIGN}%s, %f\N{DEGREE SIGN}%s' 

170 % (abs(lat), ns, abs(lon), ew)) 

171 

172 def set_longitude_grid(self, degrees): 

173 """ 

174 Set the number of degrees between each longitude grid. 

175 """ 

176 # Skip -180 and 180, which are the fixed limits. 

177 grid = np.arange(-180 + degrees, 180, degrees) 

178 self.xaxis.set_major_locator(FixedLocator(np.deg2rad(grid))) 

179 self.xaxis.set_major_formatter(self.ThetaFormatter(degrees)) 

180 

181 def set_latitude_grid(self, degrees): 

182 """ 

183 Set the number of degrees between each latitude grid. 

184 """ 

185 # Skip -90 and 90, which are the fixed limits. 

186 grid = np.arange(-90 + degrees, 90, degrees) 

187 self.yaxis.set_major_locator(FixedLocator(np.deg2rad(grid))) 

188 self.yaxis.set_major_formatter(self.ThetaFormatter(degrees)) 

189 

190 def set_longitude_grid_ends(self, degrees): 

191 """ 

192 Set the latitude(s) at which to stop drawing the longitude grids. 

193 """ 

194 self._longitude_cap = np.deg2rad(degrees) 

195 self._xaxis_pretransform \ 

196 .clear() \ 

197 .scale(1.0, self._longitude_cap * 2.0) \ 

198 .translate(0.0, -self._longitude_cap) 

199 

200 def get_data_ratio(self): 

201 ''' 

202 Return the aspect ratio of the data itself. 

203 ''' 

204 return 1.0 

205 

206 ### Interactive panning 

207 

208 def can_zoom(self): 

209 """ 

210 Return *True* if this axes supports the zoom box button functionality. 

211 

212 This axes object does not support interactive zoom box. 

213 """ 

214 return False 

215 

216 def can_pan(self): 

217 """ 

218 Return *True* if this axes supports the pan/zoom button functionality. 

219 

220 This axes object does not support interactive pan/zoom. 

221 """ 

222 return False 

223 

224 def start_pan(self, x, y, button): 

225 pass 

226 

227 def end_pan(self): 

228 pass 

229 

230 def drag_pan(self, button, key, x, y): 

231 pass 

232 

233 

234class _GeoTransform(Transform): 

235 # Factoring out some common functionality. 

236 input_dims = 2 

237 output_dims = 2 

238 is_separable = False 

239 

240 def __init__(self, resolution): 

241 """ 

242 Create a new geographical transform. 

243 

244 Resolution is the number of steps to interpolate between each input 

245 line segment to approximate its path in curved space. 

246 """ 

247 Transform.__init__(self) 

248 self._resolution = resolution 

249 

250 def __str__(self): 

251 return "{}({})".format(type(self).__name__, self._resolution) 

252 

253 def transform_path_non_affine(self, path): 

254 # docstring inherited 

255 ipath = path.interpolated(self._resolution) 

256 return Path(self.transform(ipath.vertices), ipath.codes) 

257 

258 

259class AitoffAxes(GeoAxes): 

260 name = 'aitoff' 

261 

262 class AitoffTransform(_GeoTransform): 

263 """The base Aitoff transform.""" 

264 

265 def transform_non_affine(self, ll): 

266 # docstring inherited 

267 longitude, latitude = ll.T 

268 

269 # Pre-compute some values 

270 half_long = longitude / 2.0 

271 cos_latitude = np.cos(latitude) 

272 

273 alpha = np.arccos(cos_latitude * np.cos(half_long)) 

274 # Avoid divide-by-zero errors using same method as NumPy. 

275 alpha[alpha == 0.0] = 1e-20 

276 # We want unnormalized sinc. numpy.sinc gives us normalized 

277 sinc_alpha = np.sin(alpha) / alpha 

278 

279 x = (cos_latitude * np.sin(half_long)) / sinc_alpha 

280 y = np.sin(latitude) / sinc_alpha 

281 return np.column_stack([x, y]) 

282 

283 def inverted(self): 

284 # docstring inherited 

285 return AitoffAxes.InvertedAitoffTransform(self._resolution) 

286 

287 class InvertedAitoffTransform(_GeoTransform): 

288 

289 def transform_non_affine(self, xy): 

290 # docstring inherited 

291 # MGDTODO: Math is hard ;( 

292 return xy 

293 

294 def inverted(self): 

295 # docstring inherited 

296 return AitoffAxes.AitoffTransform(self._resolution) 

297 

298 def __init__(self, *args, **kwargs): 

299 self._longitude_cap = np.pi / 2.0 

300 GeoAxes.__init__(self, *args, **kwargs) 

301 self.set_aspect(0.5, adjustable='box', anchor='C') 

302 self.cla() 

303 

304 def _get_core_transform(self, resolution): 

305 return self.AitoffTransform(resolution) 

306 

307 

308class HammerAxes(GeoAxes): 

309 name = 'hammer' 

310 

311 class HammerTransform(_GeoTransform): 

312 """The base Hammer transform.""" 

313 

314 def transform_non_affine(self, ll): 

315 # docstring inherited 

316 longitude, latitude = ll.T 

317 half_long = longitude / 2.0 

318 cos_latitude = np.cos(latitude) 

319 sqrt2 = np.sqrt(2.0) 

320 alpha = np.sqrt(1.0 + cos_latitude * np.cos(half_long)) 

321 x = (2.0 * sqrt2) * (cos_latitude * np.sin(half_long)) / alpha 

322 y = (sqrt2 * np.sin(latitude)) / alpha 

323 return np.column_stack([x, y]) 

324 

325 def inverted(self): 

326 # docstring inherited 

327 return HammerAxes.InvertedHammerTransform(self._resolution) 

328 

329 class InvertedHammerTransform(_GeoTransform): 

330 

331 def transform_non_affine(self, xy): 

332 # docstring inherited 

333 x, y = xy.T 

334 z = np.sqrt(1 - (x / 4) ** 2 - (y / 2) ** 2) 

335 longitude = 2 * np.arctan((z * x) / (2 * (2 * z ** 2 - 1))) 

336 latitude = np.arcsin(y*z) 

337 return np.column_stack([longitude, latitude]) 

338 

339 def inverted(self): 

340 # docstring inherited 

341 return HammerAxes.HammerTransform(self._resolution) 

342 

343 def __init__(self, *args, **kwargs): 

344 self._longitude_cap = np.pi / 2.0 

345 GeoAxes.__init__(self, *args, **kwargs) 

346 self.set_aspect(0.5, adjustable='box', anchor='C') 

347 self.cla() 

348 

349 def _get_core_transform(self, resolution): 

350 return self.HammerTransform(resolution) 

351 

352 

353class MollweideAxes(GeoAxes): 

354 name = 'mollweide' 

355 

356 class MollweideTransform(_GeoTransform): 

357 """The base Mollweide transform.""" 

358 

359 def transform_non_affine(self, ll): 

360 # docstring inherited 

361 def d(theta): 

362 delta = (-(theta + np.sin(theta) - pi_sin_l) 

363 / (1 + np.cos(theta))) 

364 return delta, np.abs(delta) > 0.001 

365 

366 longitude, latitude = ll.T 

367 

368 clat = np.pi/2 - np.abs(latitude) 

369 ihigh = clat < 0.087 # within 5 degrees of the poles 

370 ilow = ~ihigh 

371 aux = np.empty(latitude.shape, dtype=float) 

372 

373 if ilow.any(): # Newton-Raphson iteration 

374 pi_sin_l = np.pi * np.sin(latitude[ilow]) 

375 theta = 2.0 * latitude[ilow] 

376 delta, large_delta = d(theta) 

377 while np.any(large_delta): 

378 theta[large_delta] += delta[large_delta] 

379 delta, large_delta = d(theta) 

380 aux[ilow] = theta / 2 

381 

382 if ihigh.any(): # Taylor series-based approx. solution 

383 e = clat[ihigh] 

384 d = 0.5 * (3 * np.pi * e**2) ** (1.0/3) 

385 aux[ihigh] = (np.pi/2 - d) * np.sign(latitude[ihigh]) 

386 

387 xy = np.empty(ll.shape, dtype=float) 

388 xy[:, 0] = (2.0 * np.sqrt(2.0) / np.pi) * longitude * np.cos(aux) 

389 xy[:, 1] = np.sqrt(2.0) * np.sin(aux) 

390 

391 return xy 

392 

393 def inverted(self): 

394 # docstring inherited 

395 return MollweideAxes.InvertedMollweideTransform(self._resolution) 

396 

397 class InvertedMollweideTransform(_GeoTransform): 

398 

399 def transform_non_affine(self, xy): 

400 # docstring inherited 

401 x, y = xy.T 

402 # from Equations (7, 8) of 

403 # http://mathworld.wolfram.com/MollweideProjection.html 

404 theta = np.arcsin(y / np.sqrt(2)) 

405 longitude = (np.pi / (2 * np.sqrt(2))) * x / np.cos(theta) 

406 latitude = np.arcsin((2 * theta + np.sin(2 * theta)) / np.pi) 

407 return np.column_stack([longitude, latitude]) 

408 

409 def inverted(self): 

410 # docstring inherited 

411 return MollweideAxes.MollweideTransform(self._resolution) 

412 

413 def __init__(self, *args, **kwargs): 

414 self._longitude_cap = np.pi / 2.0 

415 GeoAxes.__init__(self, *args, **kwargs) 

416 self.set_aspect(0.5, adjustable='box', anchor='C') 

417 self.cla() 

418 

419 def _get_core_transform(self, resolution): 

420 return self.MollweideTransform(resolution) 

421 

422 

423class LambertAxes(GeoAxes): 

424 name = 'lambert' 

425 

426 class LambertTransform(_GeoTransform): 

427 """The base Lambert transform.""" 

428 

429 def __init__(self, center_longitude, center_latitude, resolution): 

430 """ 

431 Create a new Lambert transform. Resolution is the number of steps 

432 to interpolate between each input line segment to approximate its 

433 path in curved Lambert space. 

434 """ 

435 _GeoTransform.__init__(self, resolution) 

436 self._center_longitude = center_longitude 

437 self._center_latitude = center_latitude 

438 

439 def transform_non_affine(self, ll): 

440 # docstring inherited 

441 longitude, latitude = ll.T 

442 clong = self._center_longitude 

443 clat = self._center_latitude 

444 cos_lat = np.cos(latitude) 

445 sin_lat = np.sin(latitude) 

446 diff_long = longitude - clong 

447 cos_diff_long = np.cos(diff_long) 

448 

449 inner_k = np.maximum( # Prevent divide-by-zero problems 

450 1 + np.sin(clat)*sin_lat + np.cos(clat)*cos_lat*cos_diff_long, 

451 1e-15) 

452 k = np.sqrt(2 / inner_k) 

453 x = k * cos_lat*np.sin(diff_long) 

454 y = k * (np.cos(clat)*sin_lat - np.sin(clat)*cos_lat*cos_diff_long) 

455 

456 return np.column_stack([x, y]) 

457 

458 def inverted(self): 

459 # docstring inherited 

460 return LambertAxes.InvertedLambertTransform( 

461 self._center_longitude, 

462 self._center_latitude, 

463 self._resolution) 

464 

465 class InvertedLambertTransform(_GeoTransform): 

466 

467 def __init__(self, center_longitude, center_latitude, resolution): 

468 _GeoTransform.__init__(self, resolution) 

469 self._center_longitude = center_longitude 

470 self._center_latitude = center_latitude 

471 

472 def transform_non_affine(self, xy): 

473 # docstring inherited 

474 x, y = xy.T 

475 clong = self._center_longitude 

476 clat = self._center_latitude 

477 p = np.maximum(np.hypot(x, y), 1e-9) 

478 c = 2 * np.arcsin(0.5 * p) 

479 sin_c = np.sin(c) 

480 cos_c = np.cos(c) 

481 

482 latitude = np.arcsin(cos_c*np.sin(clat) + 

483 ((y*sin_c*np.cos(clat)) / p)) 

484 longitude = clong + np.arctan( 

485 (x*sin_c) / (p*np.cos(clat)*cos_c - y*np.sin(clat)*sin_c)) 

486 

487 return np.column_stack([longitude, latitude]) 

488 

489 def inverted(self): 

490 # docstring inherited 

491 return LambertAxes.LambertTransform( 

492 self._center_longitude, 

493 self._center_latitude, 

494 self._resolution) 

495 

496 def __init__(self, *args, center_longitude=0, center_latitude=0, **kwargs): 

497 self._longitude_cap = np.pi / 2 

498 self._center_longitude = center_longitude 

499 self._center_latitude = center_latitude 

500 GeoAxes.__init__(self, *args, **kwargs) 

501 self.set_aspect('equal', adjustable='box', anchor='C') 

502 self.cla() 

503 

504 def cla(self): 

505 GeoAxes.cla(self) 

506 self.yaxis.set_major_formatter(NullFormatter()) 

507 

508 def _get_core_transform(self, resolution): 

509 return self.LambertTransform( 

510 self._center_longitude, 

511 self._center_latitude, 

512 resolution) 

513 

514 def _get_affine_transform(self): 

515 return Affine2D() \ 

516 .scale(0.25) \ 

517 .translate(0.5, 0.5)