Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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, clear=False)
34 self.yaxis = maxis.YAxis(self, clear=False)
35 self.spines['geo'].register_axis(self.yaxis)
37 def clear(self):
38 # docstring inherited
39 super().clear()
41 self.set_longitude_grid(30)
42 self.set_latitude_grid(15)
43 self.set_longitude_grid_ends(75)
44 self.xaxis.set_minor_locator(NullLocator())
45 self.yaxis.set_minor_locator(NullLocator())
46 self.xaxis.set_ticks_position('none')
47 self.yaxis.set_ticks_position('none')
48 self.yaxis.set_tick_params(label1On=True)
49 # Why do we need to turn on yaxis tick labels, but
50 # xaxis tick labels are already on?
52 self.grid(mpl.rcParams['axes.grid'])
54 Axes.set_xlim(self, -np.pi, np.pi)
55 Axes.set_ylim(self, -np.pi / 2.0, np.pi / 2.0)
57 def _set_lim_and_transforms(self):
58 # A (possibly non-linear) projection on the (already scaled) data
59 self.transProjection = self._get_core_transform(self.RESOLUTION)
61 self.transAffine = self._get_affine_transform()
63 self.transAxes = BboxTransformTo(self.bbox)
65 # The complete data transformation stack -- from data all the
66 # way to display coordinates
67 self.transData = \
68 self.transProjection + \
69 self.transAffine + \
70 self.transAxes
72 # This is the transform for longitude ticks.
73 self._xaxis_pretransform = \
74 Affine2D() \
75 .scale(1, self._longitude_cap * 2) \
76 .translate(0, -self._longitude_cap)
77 self._xaxis_transform = \
78 self._xaxis_pretransform + \
79 self.transData
80 self._xaxis_text1_transform = \
81 Affine2D().scale(1, 0) + \
82 self.transData + \
83 Affine2D().translate(0, 4)
84 self._xaxis_text2_transform = \
85 Affine2D().scale(1, 0) + \
86 self.transData + \
87 Affine2D().translate(0, -4)
89 # This is the transform for latitude ticks.
90 yaxis_stretch = Affine2D().scale(np.pi * 2, 1).translate(-np.pi, 0)
91 yaxis_space = Affine2D().scale(1, 1.1)
92 self._yaxis_transform = \
93 yaxis_stretch + \
94 self.transData
95 yaxis_text_base = \
96 yaxis_stretch + \
97 self.transProjection + \
98 (yaxis_space +
99 self.transAffine +
100 self.transAxes)
101 self._yaxis_text1_transform = \
102 yaxis_text_base + \
103 Affine2D().translate(-8, 0)
104 self._yaxis_text2_transform = \
105 yaxis_text_base + \
106 Affine2D().translate(8, 0)
108 def _get_affine_transform(self):
109 transform = self._get_core_transform(1)
110 xscale, _ = transform.transform((np.pi, 0))
111 _, yscale = transform.transform((0, np.pi/2))
112 return Affine2D() \
113 .scale(0.5 / xscale, 0.5 / yscale) \
114 .translate(0.5, 0.5)
116 def get_xaxis_transform(self, which='grid'):
117 _api.check_in_list(['tick1', 'tick2', 'grid'], which=which)
118 return self._xaxis_transform
120 def get_xaxis_text1_transform(self, pad):
121 return self._xaxis_text1_transform, 'bottom', 'center'
123 def get_xaxis_text2_transform(self, pad):
124 return self._xaxis_text2_transform, 'top', 'center'
126 def get_yaxis_transform(self, which='grid'):
127 _api.check_in_list(['tick1', 'tick2', 'grid'], which=which)
128 return self._yaxis_transform
130 def get_yaxis_text1_transform(self, pad):
131 return self._yaxis_text1_transform, 'center', 'right'
133 def get_yaxis_text2_transform(self, pad):
134 return self._yaxis_text2_transform, 'center', 'left'
136 def _gen_axes_patch(self):
137 return Circle((0.5, 0.5), 0.5)
139 def _gen_axes_spines(self):
140 return {'geo': mspines.Spine.circular_spine(self, (0.5, 0.5), 0.5)}
142 def set_yscale(self, *args, **kwargs):
143 if args[0] != 'linear':
144 raise NotImplementedError
146 set_xscale = set_yscale
148 def set_xlim(self, *args, **kwargs):
149 """Not supported. Please consider using Cartopy."""
150 raise TypeError("Changing axes limits of a geographic projection is "
151 "not supported. Please consider using Cartopy.")
153 set_ylim = set_xlim
155 def format_coord(self, lon, lat):
156 """Return a format string formatting the coordinate."""
157 lon, lat = np.rad2deg([lon, lat])
158 ns = 'N' if lat >= 0.0 else 'S'
159 ew = 'E' if lon >= 0.0 else 'W'
160 return ('%f\N{DEGREE SIGN}%s, %f\N{DEGREE SIGN}%s'
161 % (abs(lat), ns, abs(lon), ew))
163 def set_longitude_grid(self, degrees):
164 """
165 Set the number of degrees between each longitude grid.
166 """
167 # Skip -180 and 180, which are the fixed limits.
168 grid = np.arange(-180 + degrees, 180, degrees)
169 self.xaxis.set_major_locator(FixedLocator(np.deg2rad(grid)))
170 self.xaxis.set_major_formatter(self.ThetaFormatter(degrees))
172 def set_latitude_grid(self, degrees):
173 """
174 Set the number of degrees between each latitude grid.
175 """
176 # Skip -90 and 90, which are the fixed limits.
177 grid = np.arange(-90 + degrees, 90, degrees)
178 self.yaxis.set_major_locator(FixedLocator(np.deg2rad(grid)))
179 self.yaxis.set_major_formatter(self.ThetaFormatter(degrees))
181 def set_longitude_grid_ends(self, degrees):
182 """
183 Set the latitude(s) at which to stop drawing the longitude grids.
184 """
185 self._longitude_cap = np.deg2rad(degrees)
186 self._xaxis_pretransform \
187 .clear() \
188 .scale(1.0, self._longitude_cap * 2.0) \
189 .translate(0.0, -self._longitude_cap)
191 def get_data_ratio(self):
192 """Return the aspect ratio of the data itself."""
193 return 1.0
195 ### Interactive panning
197 def can_zoom(self):
198 """
199 Return whether this Axes supports the zoom box button functionality.
201 This Axes object does not support interactive zoom box.
202 """
203 return False
205 def can_pan(self):
206 """
207 Return whether this Axes supports the pan/zoom button functionality.
209 This Axes object does not support interactive pan/zoom.
210 """
211 return False
213 def start_pan(self, x, y, button):
214 pass
216 def end_pan(self):
217 pass
219 def drag_pan(self, button, key, x, y):
220 pass
223class _GeoTransform(Transform):
224 # Factoring out some common functionality.
225 input_dims = output_dims = 2
227 def __init__(self, resolution):
228 """
229 Create a new geographical transform.
231 Resolution is the number of steps to interpolate between each input
232 line segment to approximate its path in curved space.
233 """
234 super().__init__()
235 self._resolution = resolution
237 def __str__(self):
238 return f"{type(self).__name__}({self._resolution})"
240 def transform_path_non_affine(self, path):
241 # docstring inherited
242 ipath = path.interpolated(self._resolution)
243 return Path(self.transform(ipath.vertices), ipath.codes)
246class AitoffAxes(GeoAxes):
247 name = 'aitoff'
249 class AitoffTransform(_GeoTransform):
250 """The base Aitoff transform."""
252 @_api.rename_parameter("3.8", "ll", "values")
253 def transform_non_affine(self, values):
254 # docstring inherited
255 longitude, latitude = values.T
257 # Pre-compute some values
258 half_long = longitude / 2.0
259 cos_latitude = np.cos(latitude)
261 alpha = np.arccos(cos_latitude * np.cos(half_long))
262 sinc_alpha = np.sinc(alpha / np.pi) # np.sinc is sin(pi*x)/(pi*x).
264 x = (cos_latitude * np.sin(half_long)) / sinc_alpha
265 y = np.sin(latitude) / sinc_alpha
266 return np.column_stack([x, y])
268 def inverted(self):
269 # docstring inherited
270 return AitoffAxes.InvertedAitoffTransform(self._resolution)
272 class InvertedAitoffTransform(_GeoTransform):
274 @_api.rename_parameter("3.8", "xy", "values")
275 def transform_non_affine(self, values):
276 # docstring inherited
277 # MGDTODO: Math is hard ;(
278 return np.full_like(values, np.nan)
280 def inverted(self):
281 # docstring inherited
282 return AitoffAxes.AitoffTransform(self._resolution)
284 def __init__(self, *args, **kwargs):
285 self._longitude_cap = np.pi / 2.0
286 super().__init__(*args, **kwargs)
287 self.set_aspect(0.5, adjustable='box', anchor='C')
288 self.clear()
290 def _get_core_transform(self, resolution):
291 return self.AitoffTransform(resolution)
294class HammerAxes(GeoAxes):
295 name = 'hammer'
297 class HammerTransform(_GeoTransform):
298 """The base Hammer transform."""
300 @_api.rename_parameter("3.8", "ll", "values")
301 def transform_non_affine(self, values):
302 # docstring inherited
303 longitude, latitude = values.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])
312 def inverted(self):
313 # docstring inherited
314 return HammerAxes.InvertedHammerTransform(self._resolution)
316 class InvertedHammerTransform(_GeoTransform):
318 @_api.rename_parameter("3.8", "xy", "values")
319 def transform_non_affine(self, values):
320 # docstring inherited
321 x, y = values.T
322 z = np.sqrt(1 - (x / 4) ** 2 - (y / 2) ** 2)
323 longitude = 2 * np.arctan((z * x) / (2 * (2 * z ** 2 - 1)))
324 latitude = np.arcsin(y*z)
325 return np.column_stack([longitude, latitude])
327 def inverted(self):
328 # docstring inherited
329 return HammerAxes.HammerTransform(self._resolution)
331 def __init__(self, *args, **kwargs):
332 self._longitude_cap = np.pi / 2.0
333 super().__init__(*args, **kwargs)
334 self.set_aspect(0.5, adjustable='box', anchor='C')
335 self.clear()
337 def _get_core_transform(self, resolution):
338 return self.HammerTransform(resolution)
341class MollweideAxes(GeoAxes):
342 name = 'mollweide'
344 class MollweideTransform(_GeoTransform):
345 """The base Mollweide transform."""
347 @_api.rename_parameter("3.8", "ll", "values")
348 def transform_non_affine(self, values):
349 # docstring inherited
350 def d(theta):
351 delta = (-(theta + np.sin(theta) - pi_sin_l)
352 / (1 + np.cos(theta)))
353 return delta, np.abs(delta) > 0.001
355 longitude, latitude = values.T
357 clat = np.pi/2 - np.abs(latitude)
358 ihigh = clat < 0.087 # within 5 degrees of the poles
359 ilow = ~ihigh
360 aux = np.empty(latitude.shape, dtype=float)
362 if ilow.any(): # Newton-Raphson iteration
363 pi_sin_l = np.pi * np.sin(latitude[ilow])
364 theta = 2.0 * latitude[ilow]
365 delta, large_delta = d(theta)
366 while np.any(large_delta):
367 theta[large_delta] += delta[large_delta]
368 delta, large_delta = d(theta)
369 aux[ilow] = theta / 2
371 if ihigh.any(): # Taylor series-based approx. solution
372 e = clat[ihigh]
373 d = 0.5 * (3 * np.pi * e**2) ** (1.0/3)
374 aux[ihigh] = (np.pi/2 - d) * np.sign(latitude[ihigh])
376 xy = np.empty(values.shape, dtype=float)
377 xy[:, 0] = (2.0 * np.sqrt(2.0) / np.pi) * longitude * np.cos(aux)
378 xy[:, 1] = np.sqrt(2.0) * np.sin(aux)
380 return xy
382 def inverted(self):
383 # docstring inherited
384 return MollweideAxes.InvertedMollweideTransform(self._resolution)
386 class InvertedMollweideTransform(_GeoTransform):
388 @_api.rename_parameter("3.8", "xy", "values")
389 def transform_non_affine(self, values):
390 # docstring inherited
391 x, y = values.T
392 # from Equations (7, 8) of
393 # https://mathworld.wolfram.com/MollweideProjection.html
394 theta = np.arcsin(y / np.sqrt(2))
395 longitude = (np.pi / (2 * np.sqrt(2))) * x / np.cos(theta)
396 latitude = np.arcsin((2 * theta + np.sin(2 * theta)) / np.pi)
397 return np.column_stack([longitude, latitude])
399 def inverted(self):
400 # docstring inherited
401 return MollweideAxes.MollweideTransform(self._resolution)
403 def __init__(self, *args, **kwargs):
404 self._longitude_cap = np.pi / 2.0
405 super().__init__(*args, **kwargs)
406 self.set_aspect(0.5, adjustable='box', anchor='C')
407 self.clear()
409 def _get_core_transform(self, resolution):
410 return self.MollweideTransform(resolution)
413class LambertAxes(GeoAxes):
414 name = 'lambert'
416 class LambertTransform(_GeoTransform):
417 """The base Lambert transform."""
419 def __init__(self, center_longitude, center_latitude, resolution):
420 """
421 Create a new Lambert transform. Resolution is the number of steps
422 to interpolate between each input line segment to approximate its
423 path in curved Lambert space.
424 """
425 _GeoTransform.__init__(self, resolution)
426 self._center_longitude = center_longitude
427 self._center_latitude = center_latitude
429 @_api.rename_parameter("3.8", "ll", "values")
430 def transform_non_affine(self, values):
431 # docstring inherited
432 longitude, latitude = values.T
433 clong = self._center_longitude
434 clat = self._center_latitude
435 cos_lat = np.cos(latitude)
436 sin_lat = np.sin(latitude)
437 diff_long = longitude - clong
438 cos_diff_long = np.cos(diff_long)
440 inner_k = np.maximum( # Prevent divide-by-zero problems
441 1 + np.sin(clat)*sin_lat + np.cos(clat)*cos_lat*cos_diff_long,
442 1e-15)
443 k = np.sqrt(2 / inner_k)
444 x = k * cos_lat*np.sin(diff_long)
445 y = k * (np.cos(clat)*sin_lat - np.sin(clat)*cos_lat*cos_diff_long)
447 return np.column_stack([x, y])
449 def inverted(self):
450 # docstring inherited
451 return LambertAxes.InvertedLambertTransform(
452 self._center_longitude,
453 self._center_latitude,
454 self._resolution)
456 class InvertedLambertTransform(_GeoTransform):
458 def __init__(self, center_longitude, center_latitude, resolution):
459 _GeoTransform.__init__(self, resolution)
460 self._center_longitude = center_longitude
461 self._center_latitude = center_latitude
463 @_api.rename_parameter("3.8", "xy", "values")
464 def transform_non_affine(self, values):
465 # docstring inherited
466 x, y = values.T
467 clong = self._center_longitude
468 clat = self._center_latitude
469 p = np.maximum(np.hypot(x, y), 1e-9)
470 c = 2 * np.arcsin(0.5 * p)
471 sin_c = np.sin(c)
472 cos_c = np.cos(c)
474 latitude = np.arcsin(cos_c*np.sin(clat) +
475 ((y*sin_c*np.cos(clat)) / p))
476 longitude = clong + np.arctan(
477 (x*sin_c) / (p*np.cos(clat)*cos_c - y*np.sin(clat)*sin_c))
479 return np.column_stack([longitude, latitude])
481 def inverted(self):
482 # docstring inherited
483 return LambertAxes.LambertTransform(
484 self._center_longitude,
485 self._center_latitude,
486 self._resolution)
488 def __init__(self, *args, center_longitude=0, center_latitude=0, **kwargs):
489 self._longitude_cap = np.pi / 2
490 self._center_longitude = center_longitude
491 self._center_latitude = center_latitude
492 super().__init__(*args, **kwargs)
493 self.set_aspect('equal', adjustable='box', anchor='C')
494 self.clear()
496 def clear(self):
497 # docstring inherited
498 super().clear()
499 self.yaxis.set_major_formatter(NullFormatter())
501 def _get_core_transform(self, resolution):
502 return self.LambertTransform(
503 self._center_longitude,
504 self._center_latitude,
505 resolution)
507 def _get_affine_transform(self):
508 return Affine2D() \
509 .scale(0.25) \
510 .translate(0.5, 0.5)