Coverage for /usr/lib/python3/dist-packages/matplotlib/projections/geo.py: 33%

1import numpy as np

3import matplotlib as mpl

4from matplotlib import _api

5from matplotlib.axes import Axes

6import matplotlib.axis as maxis

7from matplotlib.patches import Circle

8from matplotlib.path import Path

9import matplotlib.spines as mspines

10from matplotlib.ticker import (

11 Formatter, NullLocator, FixedLocator, NullFormatter)

12from matplotlib.transforms import Affine2D, BboxTransformTo, Transform

15class GeoAxes(Axes):

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

18 class ThetaFormatter(Formatter):

19 """

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

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

22 """

23 def __init__(self, round_to=1.0):

24 self._round_to = round_to

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

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

28 return f"{degrees:0.0f}\N{{DEGREE SIGN}}"

30 RESOLUTION = 75

32 def _init_axis(self):

33 self.xaxis = maxis.XAxis(self)

34 self.yaxis = maxis.YAxis(self)

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

36 # Axes._init_axis() -- until GeoAxes.xaxis.clear() works.

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

38 self._update_transScale()

40 def clear(self):

41 # docstring inherited

42 super().clear()

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?

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

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

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

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)

64 self.transAffine = self._get_affine_transform()

66 self.transAxes = BboxTransformTo(self.bbox)

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

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)

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 _api.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 _api.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 """Not supported. Please consider using Cartopy."""

153 raise TypeError("Changing axes limits of a geographic projection is "

154 "not supported. Please consider using 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 ns = 'N' if lat >= 0.0 else 'S'

162 ew = 'E' if lon >= 0.0 else 'W'

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

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

165

166 def set_longitude_grid(self, degrees):

167 """

168 Set the number of degrees between each longitude grid.

169 """

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

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

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

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

174

175 def set_latitude_grid(self, degrees):

176 """

177 Set the number of degrees between each latitude grid.

178 """

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

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

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

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

183

184 def set_longitude_grid_ends(self, degrees):

185 """

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

187 """

188 self._longitude_cap = np.deg2rad(degrees)

189 self._xaxis_pretransform \

190 .clear() \

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

192 .translate(0.0, -self._longitude_cap)

193

194 def get_data_ratio(self):

195 """Return the aspect ratio of the data itself."""

196 return 1.0

197

198 ### Interactive panning

199

200 def can_zoom(self):

201 """

202 Return whether this axes supports the zoom box button functionality.

203

204 This axes object does not support interactive zoom box.

205 """

206 return False

207

208 def can_pan(self):

209 """

210 Return whether this axes supports the pan/zoom button functionality.

211

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

213 """

214 return False

215

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

217 pass

218

219 def end_pan(self):

220 pass

221

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

223 pass

224

225

226class _GeoTransform(Transform):

227 # Factoring out some common functionality.

228 input_dims = output_dims = 2

229

230 def __init__(self, resolution):

231 """

232 Create a new geographical transform.

233

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

235 line segment to approximate its path in curved space.

236 """

237 super().__init__()

238 self._resolution = resolution

239

240 def __str__(self):

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

242

243 def transform_path_non_affine(self, path):

244 # docstring inherited

245 ipath = path.interpolated(self._resolution)

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

247

248

249class AitoffAxes(GeoAxes):

250 name = 'aitoff'

251

252 class AitoffTransform(_GeoTransform):

253 """The base Aitoff transform."""

254

255 def transform_non_affine(self, ll):

256 # docstring inherited

257 longitude, latitude = ll.T

258

259 # Pre-compute some values

260 half_long = longitude / 2.0

261 cos_latitude = np.cos(latitude)

262

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

264 sinc_alpha = np.sinc(alpha / np.pi) # np.sinc is sin(pi*x)/(pi*x).

265

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

267 y = np.sin(latitude) / sinc_alpha

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

269

270 def inverted(self):

271 # docstring inherited

272 return AitoffAxes.InvertedAitoffTransform(self._resolution)

273

274 class InvertedAitoffTransform(_GeoTransform):

275

276 def transform_non_affine(self, xy):

277 # docstring inherited

278 # MGDTODO: Math is hard ;(

279 return np.full_like(xy, np.nan)

280

281 def inverted(self):

282 # docstring inherited

283 return AitoffAxes.AitoffTransform(self._resolution)

284

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

286 self._longitude_cap = np.pi / 2.0

287 super().__init__(*args, **kwargs)

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

289 self.clear()

290

291 def _get_core_transform(self, resolution):

292 return self.AitoffTransform(resolution)

293

294

295class HammerAxes(GeoAxes):

296 name = 'hammer'

297

298 class HammerTransform(_GeoTransform):

299 """The base Hammer transform."""

300

301 def transform_non_affine(self, ll):

302 # docstring inherited

303 longitude, latitude = ll.T

304 half_long = longitude / 2.0

305 cos_latitude = np.cos(latitude)

306 sqrt2 = np.sqrt(2.0)

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

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

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

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

311

312 def inverted(self):

313 # docstring inherited

314 return HammerAxes.InvertedHammerTransform(self._resolution)

315

316 class InvertedHammerTransform(_GeoTransform):

317

318 def transform_non_affine(self, xy):

319 # docstring inherited

320 x, y = xy.T

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

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

323 latitude = np.arcsin(y*z)

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

325

326 def inverted(self):

327 # docstring inherited

328 return HammerAxes.HammerTransform(self._resolution)

329

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

331 self._longitude_cap = np.pi / 2.0

332 super().__init__(*args, **kwargs)

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

334 self.clear()

335

336 def _get_core_transform(self, resolution):

337 return self.HammerTransform(resolution)

338

339

340class MollweideAxes(GeoAxes):

341 name = 'mollweide'

342

343 class MollweideTransform(_GeoTransform):

344 """The base Mollweide transform."""

345

346 def transform_non_affine(self, ll):

347 # docstring inherited

348 def d(theta):

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

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

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

352

353 longitude, latitude = ll.T

354

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

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

357 ilow = ~ihigh

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

359

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

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

362 theta = 2.0 * latitude[ilow]

363 delta, large_delta = d(theta)

364 while np.any(large_delta):

365 theta[large_delta] += delta[large_delta]

366 delta, large_delta = d(theta)

367 aux[ilow] = theta / 2

368

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

370 e = clat[ihigh]

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

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

373

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

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

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

377

378 return xy

379

380 def inverted(self):

381 # docstring inherited

382 return MollweideAxes.InvertedMollweideTransform(self._resolution)

383

384 class InvertedMollweideTransform(_GeoTransform):

385

386 def transform_non_affine(self, xy):

387 # docstring inherited

388 x, y = xy.T

389 # from Equations (7, 8) of

390 # https://mathworld.wolfram.com/MollweideProjection.html

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

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

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

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

395

396 def inverted(self):

397 # docstring inherited

398 return MollweideAxes.MollweideTransform(self._resolution)

399

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

401 self._longitude_cap = np.pi / 2.0

402 super().__init__(*args, **kwargs)

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

404 self.clear()

405

406 def _get_core_transform(self, resolution):

407 return self.MollweideTransform(resolution)

408

409

410class LambertAxes(GeoAxes):

411 name = 'lambert'

412

413 class LambertTransform(_GeoTransform):

414 """The base Lambert transform."""

415

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

417 """

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

419 to interpolate between each input line segment to approximate its

420 path in curved Lambert space.

421 """

422 _GeoTransform.__init__(self, resolution)

423 self._center_longitude = center_longitude

424 self._center_latitude = center_latitude

425

426 def transform_non_affine(self, ll):

427 # docstring inherited

428 longitude, latitude = ll.T

429 clong = self._center_longitude

430 clat = self._center_latitude

431 cos_lat = np.cos(latitude)

432 sin_lat = np.sin(latitude)

433 diff_long = longitude - clong

434 cos_diff_long = np.cos(diff_long)

435

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

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

438 1e-15)

439 k = np.sqrt(2 / inner_k)

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

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

442

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

444

445 def inverted(self):

446 # docstring inherited

447 return LambertAxes.InvertedLambertTransform(

448 self._center_longitude,

449 self._center_latitude,

450 self._resolution)

451

452 class InvertedLambertTransform(_GeoTransform):

453

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

455 _GeoTransform.__init__(self, resolution)

456 self._center_longitude = center_longitude

457 self._center_latitude = center_latitude

458

459 def transform_non_affine(self, xy):

460 # docstring inherited

461 x, y = xy.T

462 clong = self._center_longitude

463 clat = self._center_latitude

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

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

466 sin_c = np.sin(c)

467 cos_c = np.cos(c)

468

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

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

471 longitude = clong + np.arctan(

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

473

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

475

476 def inverted(self):

477 # docstring inherited

478 return LambertAxes.LambertTransform(

479 self._center_longitude,

480 self._center_latitude,

481 self._resolution)

482

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

484 self._longitude_cap = np.pi / 2

485 self._center_longitude = center_longitude

486 self._center_latitude = center_latitude

487 super().__init__(*args, **kwargs)

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

489 self.clear()

490

491 def clear(self):

492 # docstring inherited

493 super().clear()

494 self.yaxis.set_major_formatter(NullFormatter())

495

496 def _get_core_transform(self, resolution):

497 return self.LambertTransform(

498 self._center_longitude,

499 self._center_latitude,

500 resolution)

501

502 def _get_affine_transform(self):

503 return Affine2D() \

504 .scale(0.25) \

505 .translate(0.5, 0.5)