Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/matplotlib/projections/geo.py : 32%

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
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
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
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
31 RESOLUTION = 75
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()
41 def cla(self):
42 Axes.cla(self)
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(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)
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)
119 def get_xaxis_transform(self, which='grid'):
120 cbook._check_in_list(['tick1', 'tick2', 'grid'], which=which)
121 return self._xaxis_transform
123 def get_xaxis_text1_transform(self, pad):
124 return self._xaxis_text1_transform, 'bottom', 'center'
126 def get_xaxis_text2_transform(self, pad):
127 return self._xaxis_text2_transform, 'top', 'center'
129 def get_yaxis_transform(self, which='grid'):
130 cbook._check_in_list(['tick1', 'tick2', 'grid'], which=which)
131 return self._yaxis_transform
133 def get_yaxis_text1_transform(self, pad):
134 return self._yaxis_text1_transform, 'center', 'right'
136 def get_yaxis_text2_transform(self, pad):
137 return self._yaxis_text2_transform, 'center', 'left'
139 def _gen_axes_patch(self):
140 return Circle((0.5, 0.5), 0.5)
142 def _gen_axes_spines(self):
143 return {'geo': mspines.Spine.circular_spine(self, (0.5, 0.5), 0.5)}
145 def set_yscale(self, *args, **kwargs):
146 if args[0] != 'linear':
147 raise NotImplementedError
149 set_xscale = set_yscale
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.")
156 set_ylim = set_xlim
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))
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))
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))
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)
200 def get_data_ratio(self):
201 '''
202 Return the aspect ratio of the data itself.
203 '''
204 return 1.0
206 ### Interactive panning
208 def can_zoom(self):
209 """
210 Return *True* if this axes supports the zoom box button functionality.
212 This axes object does not support interactive zoom box.
213 """
214 return False
216 def can_pan(self):
217 """
218 Return *True* if this axes supports the pan/zoom button functionality.
220 This axes object does not support interactive pan/zoom.
221 """
222 return False
224 def start_pan(self, x, y, button):
225 pass
227 def end_pan(self):
228 pass
230 def drag_pan(self, button, key, x, y):
231 pass
234class _GeoTransform(Transform):
235 # Factoring out some common functionality.
236 input_dims = 2
237 output_dims = 2
238 is_separable = False
240 def __init__(self, resolution):
241 """
242 Create a new geographical transform.
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
250 def __str__(self):
251 return "{}({})".format(type(self).__name__, self._resolution)
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)
259class AitoffAxes(GeoAxes):
260 name = 'aitoff'
262 class AitoffTransform(_GeoTransform):
263 """The base Aitoff transform."""
265 def transform_non_affine(self, ll):
266 # docstring inherited
267 longitude, latitude = ll.T
269 # Pre-compute some values
270 half_long = longitude / 2.0
271 cos_latitude = np.cos(latitude)
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
279 x = (cos_latitude * np.sin(half_long)) / sinc_alpha
280 y = np.sin(latitude) / sinc_alpha
281 return np.column_stack([x, y])
283 def inverted(self):
284 # docstring inherited
285 return AitoffAxes.InvertedAitoffTransform(self._resolution)
287 class InvertedAitoffTransform(_GeoTransform):
289 def transform_non_affine(self, xy):
290 # docstring inherited
291 # MGDTODO: Math is hard ;(
292 return xy
294 def inverted(self):
295 # docstring inherited
296 return AitoffAxes.AitoffTransform(self._resolution)
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()
304 def _get_core_transform(self, resolution):
305 return self.AitoffTransform(resolution)
308class HammerAxes(GeoAxes):
309 name = 'hammer'
311 class HammerTransform(_GeoTransform):
312 """The base Hammer transform."""
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])
325 def inverted(self):
326 # docstring inherited
327 return HammerAxes.InvertedHammerTransform(self._resolution)
329 class InvertedHammerTransform(_GeoTransform):
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])
339 def inverted(self):
340 # docstring inherited
341 return HammerAxes.HammerTransform(self._resolution)
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()
349 def _get_core_transform(self, resolution):
350 return self.HammerTransform(resolution)
353class MollweideAxes(GeoAxes):
354 name = 'mollweide'
356 class MollweideTransform(_GeoTransform):
357 """The base Mollweide transform."""
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
366 longitude, latitude = ll.T
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)
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
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])
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)
391 return xy
393 def inverted(self):
394 # docstring inherited
395 return MollweideAxes.InvertedMollweideTransform(self._resolution)
397 class InvertedMollweideTransform(_GeoTransform):
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])
409 def inverted(self):
410 # docstring inherited
411 return MollweideAxes.MollweideTransform(self._resolution)
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()
419 def _get_core_transform(self, resolution):
420 return self.MollweideTransform(resolution)
423class LambertAxes(GeoAxes):
424 name = 'lambert'
426 class LambertTransform(_GeoTransform):
427 """The base Lambert transform."""
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
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)
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)
456 return np.column_stack([x, y])
458 def inverted(self):
459 # docstring inherited
460 return LambertAxes.InvertedLambertTransform(
461 self._center_longitude,
462 self._center_latitude,
463 self._resolution)
465 class InvertedLambertTransform(_GeoTransform):
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
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)
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))
487 return np.column_stack([longitude, latitude])
489 def inverted(self):
490 # docstring inherited
491 return LambertAxes.LambertTransform(
492 self._center_longitude,
493 self._center_latitude,
494 self._resolution)
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()
504 def cla(self):
505 GeoAxes.cla(self)
506 self.yaxis.set_major_formatter(NullFormatter())
508 def _get_core_transform(self, resolution):
509 return self.LambertTransform(
510 self._center_longitude,
511 self._center_latitude,
512 resolution)
514 def _get_affine_transform(self):
515 return Affine2D() \
516 .scale(0.25) \
517 .translate(0.5, 0.5)