1# -*- coding: utf-8 -*- 2""" 3Module for basemap related plotting in ObsPy. 4 5:copyright: 6 The ObsPy Development Team (devs@obspy.org) 7:license: 8 GNU Lesser General Public License, Version 3 9 (https://www.gnu.org/copyleft/lesser.html) 10""" 11from __future__ import (absolute_import, division, print_function, 12 unicode_literals) 13from future.builtins import * # NOQA @UnusedWildImport 14from future.utils import native_str 15 16import datetime 17import warnings 18 19import numpy as np 20from matplotlib.colorbar import Colorbar 21from matplotlib.dates import AutoDateFormatter, AutoDateLocator, date2num 22from matplotlib import patheffects 23from matplotlib.ticker import (FormatStrFormatter, Formatter, FuncFormatter, 24 MaxNLocator) 25 26from obspy import UTCDateTime 27from obspy.core.util import (BASEMAP_VERSION, CARTOPY_VERSION, 28 MATPLOTLIB_VERSION, PROJ4_VERSION) 29from obspy.geodetics.base import mean_longitude 30 31 32if BASEMAP_VERSION: 33 from mpl_toolkits.basemap import Basemap 34 HAS_BASEMAP = True 35 if BASEMAP_VERSION < [1, 0, 4]: 36 warnings.warn("All basemap version < 1.0.4 contain a serious bug " 37 "when rendering countries and continents. ObsPy will " 38 "still work but the maps might be wrong. Please update " 39 "your basemap installation.") 40 if PROJ4_VERSION and PROJ4_VERSION[0] == 5: 41 msg = ( 42 "basemap/pyproj with proj4 version >= 5 has a bug that results in " 43 "inverted map axes. Your maps may be wrong. Please use another " 44 "version of proj4, or use cartopy. " 45 "See https://github.com/matplotlib/basemap/issues/443") 46 warnings.warn(msg) 47 if MATPLOTLIB_VERSION == [3, 0, 1] and BASEMAP_VERSION >= [1, 1, 0]: 48 msg = ( 49 "basemap and matplotlib version 3.0.1 have compataibilty issues, " 50 "please change your matplotlib version. " 51 "See https://github.com/matplotlib/basemap/issues/435") 52 warnings.warn(msg) 53else: 54 HAS_BASEMAP = False 55 56if CARTOPY_VERSION and CARTOPY_VERSION >= [0, 12, 0]: 57 import cartopy.crs as ccrs 58 import cartopy.feature as cfeature 59 HAS_CARTOPY = True 60else: 61 HAS_CARTOPY = False 62 63 64if not HAS_BASEMAP and not HAS_CARTOPY: 65 msg = ("Neither basemap nor cartopy installed, map plots will not work.") 66 warnings.warn(msg) 67 68 69_BASEMAP_RESOLUTIONS = { 70 '110m': 'l', 71 '50m': 'i', 72 '10m': 'f', 73 'c': 'c', 74 'l': 'l', 75 'i': 'i', 76 'h': 'h', 77 'f': 'f', 78} 79 80_CARTOPY_RESOLUTIONS = { 81 'c': '110m', 82 'l': '110m', 83 'i': '50m', 84 'h': '50m', 85 'f': '10m', 86 '110m': '110m', 87 '50m': '50m', 88 '10m': '10m', 89} 90 91if HAS_CARTOPY: 92 _CARTOPY_FEATURES = { 93 '110m': (cfeature.BORDERS, cfeature.LAND, cfeature.OCEAN), 94 } 95 96 97def plot_basemap(lons, lats, size, color, labels=None, projection='global', 98 resolution='l', continent_fill_color='0.8', 99 water_fill_color='1.0', colormap=None, colorbar=None, 100 marker="o", title=None, colorbar_ticklabel_format=None, 101 show=True, fig=None, **kwargs): # @UnusedVariable 102 """ 103 Creates a basemap plot with a data point scatter plot. 104 105 :type lons: list/tuple of floats 106 :param lons: Longitudes of the data points. 107 :type lats: list/tuple of floats 108 :param lats: Latitudes of the data points. 109 :type size: float or list/tuple of floats 110 :param size: Size of the individual points in the scatter plot. 111 :type color: list/tuple of floats (or objects that can be 112 converted to floats, like e.g. 113 :class:`~obspy.core.utcdatetime.UTCDateTime`) 114 :param color: Color information of the individual data points to be 115 used in the specified color map (e.g. origin depths, 116 origin times). 117 :type labels: list/tuple of str 118 :param labels: Annotations for the individual data points. 119 :type projection: str, optional 120 :param projection: The map projection. 121 Currently supported are: 122 123 * ``"global"`` (Will plot the whole world.) 124 * ``"ortho"`` (Will center around the mean lat/long.) 125 * ``"local"`` (Will plot around local events) 126 127 Defaults to "global". 128 :type resolution: str, optional 129 :param resolution: Resolution of the boundary database to use. Will be 130 based directly to the basemap module. Possible values are: 131 132 * ``"c"`` (crude) 133 * ``"l"`` (low) 134 * ``"i"`` (intermediate) 135 * ``"h"`` (high) 136 * ``"f"`` (full) 137 138 Defaults to ``"l"``. For compatibility, you may also specify any of the 139 Cartopy resolutions defined in :func:`plot_cartopy`. 140 :type continent_fill_color: Valid matplotlib color, optional 141 :param continent_fill_color: Color of the continents. Defaults to 142 ``"0.9"`` which is a light gray. 143 :type water_fill_color: Valid matplotlib color, optional 144 :param water_fill_color: Color of all water bodies. 145 Defaults to ``"white"``. 146 :type colormap: str, any matplotlib colormap, optional 147 :param colormap: The colormap for color-coding the events as provided 148 in `color` kwarg. 149 The event with the smallest `color` property will have the 150 color of one end of the colormap and the event with the highest 151 `color` property the color of the other end with all other events 152 in between. 153 Defaults to None which will use the default matplotlib colormap. 154 :type colorbar: bool, optional 155 :param colorbar: When left `None`, a colorbar is plotted if more than one 156 object is plotted. Using `True`/`False` the colorbar can be forced 157 on/off. 158 :type title: str 159 :param title: Title above plot. 160 :type colorbar_ticklabel_format: str or function or 161 subclass of :class:`matplotlib.ticker.Formatter` 162 :param colorbar_ticklabel_format: Format string or Formatter used to format 163 colorbar tick labels. 164 :type show: bool 165 :param show: Whether to show the figure after plotting or not. Can be used 166 to do further customization of the plot before showing it. 167 :type fig: :class:`matplotlib.figure.Figure` 168 :param fig: Figure instance to reuse, returned from a previous 169 :func:`plot_basemap` call. If a previous basemap plot is reused, any 170 kwargs regarding the basemap plot setup will be ignored (i.e. 171 `projection`, `resolution`, `continent_fill_color`, 172 `water_fill_color`). Note that multiple plots using colorbars likely 173 are problematic, but e.g. one station plot (without colorbar) and one 174 event plot (with colorbar) together should work well. 175 """ 176 import matplotlib.pyplot as plt 177 178 if any([isinstance(c, (datetime.datetime, UTCDateTime)) for c in color]): 179 datetimeplot = True 180 color = [ 181 (np.isfinite(float(t)) and 182 date2num(getattr(t, 'datetime', t)) or 183 np.nan) 184 for t in color] 185 else: 186 datetimeplot = False 187 188 # The colorbar should only be plotted if more then one event is 189 # present. 190 if colorbar is None: 191 if len(lons) > 1 and hasattr(color, "__len__") and \ 192 not isinstance(color, (str, native_str)): 193 colorbar = True 194 else: 195 colorbar = False 196 197 if fig is None: 198 fig = plt.figure() 199 200 if projection == "local": 201 ax_x0, ax_width = 0.10, 0.80 202 elif projection == "global": 203 ax_x0, ax_width = 0.01, 0.98 204 else: 205 ax_x0, ax_width = 0.05, 0.90 206 207 if colorbar: 208 map_ax = fig.add_axes([ax_x0, 0.13, ax_width, 0.77]) 209 cm_ax = fig.add_axes([ax_x0, 0.05, ax_width, 0.05]) 210 else: 211 ax_y0, ax_height = 0.05, 0.85 212 if projection == "local": 213 ax_y0 += 0.05 214 ax_height -= 0.05 215 map_ax = fig.add_axes([ax_x0, ax_y0, ax_width, ax_height]) 216 bmap = None 217 218 else: 219 error_message_suffix = ( 220 ". Please provide a figure object from a previous call to the " 221 ".plot() method of e.g. an Inventory or Catalog object.") 222 try: 223 map_ax = fig.axes[0] 224 except IndexError as e: 225 e.args = tuple([e.args[0] + error_message_suffix] + 226 list(e.args[1:])) 227 raise 228 try: 229 bmap = fig.bmap 230 except AttributeError as e: 231 e.args = tuple([e.args[0] + error_message_suffix] + 232 list(e.args[1:])) 233 raise 234 235 # basemap plots will break with basemap 1.1.0 together with matplotlib 236 # >=2.3 (see matplotlib/basemap#382) so only thing we can do is show a 237 # nicer message. 238 # XXX can be removed maybe a year or so after basemap 239 # 1.1.1 or 1.2.0 is released 240 try: 241 from matplotlib.cbook import MatplotlibDeprecationWarning 242 except ImportError: 243 # matplotlib 1.2.0 does not have that warning class yet 244 # XXX can be removed when minimum matplotlib version gets bumped to 245 # XXX 1.3.0 246 category = {} 247 else: 248 category = {'category': MatplotlibDeprecationWarning} 249 try: 250 with warnings.catch_warnings(): 251 warnings.filterwarnings( 252 'ignore', message='The axesPatch function was deprecated ' 253 'in version 2.1. Use Axes.patch instead.', 254 module='.*basemap.*', **category) 255 scatter = _plot_basemap_into_axes( 256 ax=map_ax, lons=lons, lats=lats, size=size, color=color, 257 bmap=bmap, labels=labels, projection=projection, 258 resolution=resolution, 259 continent_fill_color=continent_fill_color, 260 water_fill_color=water_fill_color, colormap=colormap, 261 marker=marker, title="", adjust_aspect_to_colorbar=colorbar, 262 **kwargs) 263 except AttributeError as e: 264 if 'axesPatch' not in str(e): 265 raise 266 msg = ('Encountered a problem doing the basemap plot due to a known ' 267 'issue of matplotlib >=2.3 together with basemap <=1.1.0 (see ' 268 'https://github.com/matplotlib/basemap/issues/382). Please ' 269 'update basemap to a version >1.1.0 if available or downgrade ' 270 'matplotlib to a version <2.3.') 271 raise Exception(msg) 272 273 if title: 274 plt.suptitle(title) 275 276 if colorbar: 277 if colorbar_ticklabel_format is not None: 278 if isinstance(colorbar_ticklabel_format, (str, native_str)): 279 formatter = FormatStrFormatter(colorbar_ticklabel_format) 280 elif hasattr(colorbar_ticklabel_format, '__call__'): 281 formatter = FuncFormatter(colorbar_ticklabel_format) 282 elif isinstance(colorbar_ticklabel_format, Formatter): 283 formatter = colorbar_ticklabel_format 284 locator = MaxNLocator(5) 285 else: 286 if datetimeplot: 287 locator = AutoDateLocator() 288 formatter = AutoDateFormatter(locator) 289 # Compat with old matplotlib versions. 290 if hasattr(formatter, "scaled"): 291 formatter.scaled[1 / (24. * 60.)] = '%H:%M:%S' 292 else: 293 locator = None 294 formatter = None 295 296 # normal case: axes for colorbar was set up in this routine 297 if "cm_ax" in locals(): 298 cb_kwargs = {"cax": cm_ax} 299 # unusual case: reusing a plot that has no colorbar set up previously 300 else: 301 cb_kwargs = {"ax": map_ax} 302 303 cb = fig.colorbar( 304 mappable=scatter, cmap=colormap, orientation='horizontal', 305 ticks=locator, format=formatter, **cb_kwargs) 306 # Compat with old matplotlib versions. 307 if hasattr(cb, "update_ticks"): 308 cb.update_ticks() 309 310 if show: 311 plt.show() 312 313 return fig 314 315 316def _plot_basemap_into_axes( 317 ax, lons, lats, size, color, bmap=None, labels=None, 318 projection='global', resolution='l', continent_fill_color='0.8', 319 water_fill_color='1.0', colormap=None, marker="o", title=None, 320 adjust_aspect_to_colorbar=False, **kwargs): # @UnusedVariable 321 """ 322 Creates a (or adds to existing) basemap plot with a data point scatter 323 plot in given axes. 324 325 See :func:`plot_basemap` for details on most args/kwargs. 326 327 :type ax: :class:`matplotlib.axes.Axes` 328 :param ax: Existing matplotlib axes instance, optionally with previous 329 basemap plot (see `bmap` kwarg). 330 :type bmap: :class:`mpl_toolkits.basemap.Basemap` 331 :param bmap: Basemap instance in provided matplotlib Axes `ax` to reuse. If 332 specified, any kwargs regarding the basemap plot setup will be ignored 333 (i.e. `projection`, `resolution`, `continent_fill_color`, 334 `water_fill_color`). 335 :rtype: :class:`matplotlib.collections.PathCollection` 336 :returns: Matplotlib path collection (e.g. to reuse for colorbars). 337 """ 338 fig = ax.figure 339 if bmap is None: 340 341 if projection == 'global': 342 bmap = Basemap(projection='moll', lon_0=round(np.mean(lons), 4), 343 resolution=_BASEMAP_RESOLUTIONS[resolution], 344 ax=ax) 345 elif projection == 'ortho': 346 bmap = Basemap(projection='ortho', 347 resolution=_BASEMAP_RESOLUTIONS[resolution], 348 lat_0=round(np.mean(lats), 4), 349 lon_0=round(mean_longitude(lons), 4), ax=ax) 350 elif projection == 'local': 351 if min(lons) < -150 and max(lons) > 150: 352 max_lons = max(np.array(lons) % 360) 353 min_lons = min(np.array(lons) % 360) 354 else: 355 max_lons = max(lons) 356 min_lons = min(lons) 357 lat_0 = max(lats) / 2. + min(lats) / 2. 358 lon_0 = max_lons / 2. + min_lons / 2. 359 if lon_0 > 180: 360 lon_0 -= 360 361 deg2m_lat = 2 * np.pi * 6371 * 1000 / 360 362 deg2m_lon = deg2m_lat * np.cos(lat_0 / 180 * np.pi) 363 if len(lats) > 1: 364 height = (max(lats) - min(lats)) * deg2m_lat 365 width = (max_lons - min_lons) * deg2m_lon 366 margin = 0.2 * (width + height) 367 height += margin 368 width += margin 369 else: 370 height = 2.0 * deg2m_lat 371 width = 5.0 * deg2m_lon 372 # do intelligent aspect calculation for local projection 373 # adjust to figure dimensions 374 w, h = fig.get_size_inches() 375 ax_bbox = ax.get_position() 376 aspect = (w * ax_bbox.width) / (h * ax_bbox.height) 377 if adjust_aspect_to_colorbar: 378 aspect *= 1.2 379 if width / height < aspect: 380 width = height * aspect 381 else: 382 height = width / aspect 383 384 bmap = Basemap(projection='aea', 385 resolution=_BASEMAP_RESOLUTIONS[resolution], 386 lat_0=round(lat_0, 4), 387 lon_0=round(lon_0, 4), 388 width=width, height=height, ax=ax) 389 # not most elegant way to calculate some round lats/lons 390 391 def linspace2(val1, val2, n): 392 """ 393 returns around n 'nice' values between val1 and val2 394 """ 395 dval = val2 - val1 396 round_pos = int(round(-np.log10(1. * dval / n))) 397 # Fake negative rounding as not supported by future as of now. 398 if round_pos < 0: 399 factor = 10 ** (abs(round_pos)) 400 delta = round(2. * dval / n / factor) * factor / 2 401 else: 402 delta = round(2. * dval / n, round_pos) / 2 403 new_val1 = np.ceil(val1 / delta) * delta 404 new_val2 = np.floor(val2 / delta) * delta 405 n = int((new_val2 - new_val1) / delta + 1) 406 return np.linspace(new_val1, new_val2, n) 407 408 n_1 = int(np.ceil(height / max(width, height) * 8)) 409 n_2 = int(np.ceil(width / max(width, height) * 8)) 410 parallels = linspace2(lat_0 - height / 2 / deg2m_lat, 411 lat_0 + height / 2 / deg2m_lat, n_1) 412 413 # Old basemap versions have problems with non-integer parallels. 414 try: 415 bmap.drawparallels(parallels, labels=[0, 1, 1, 0]) 416 except KeyError: 417 parallels = sorted(list(set(map(int, parallels)))) 418 bmap.drawparallels(parallels, labels=[0, 1, 1, 0]) 419 420 if min(lons) < -150 and max(lons) > 150: 421 lon_0 %= 360 422 meridians = linspace2(lon_0 - width / 2 / deg2m_lon, 423 lon_0 + width / 2 / deg2m_lon, n_2) 424 meridians[meridians > 180] -= 360 425 bmap.drawmeridians(meridians, labels=[1, 0, 0, 1]) 426 else: 427 msg = "Projection '%s' not supported." % projection 428 raise ValueError(msg) 429 430 # draw coast lines, country boundaries, fill continents. 431 if MATPLOTLIB_VERSION >= [2, 0, 0]: 432 ax.set_facecolor(water_fill_color) 433 else: 434 ax.set_axis_bgcolor(water_fill_color) 435 # newer matplotlib errors out if called with empty coastline data (no 436 # coast on map) 437 if np.size(getattr(bmap, 'coastsegs', [])): 438 bmap.drawcoastlines(color="0.4") 439 bmap.drawcountries(color="0.75") 440 bmap.fillcontinents(color=continent_fill_color, 441 lake_color=water_fill_color) 442 # draw the edge of the bmap projection region (the projection limb) 443 bmap.drawmapboundary(fill_color=water_fill_color) 444 # draw lat/lon grid lines every 30 degrees. 445 bmap.drawmeridians(np.arange(-180, 180, 30)) 446 bmap.drawparallels(np.arange(-90, 90, 30)) 447 fig.bmap = bmap 448 449 # compute the native bmap projection coordinates for events. 450 x, y = bmap(lons, lats) 451 # plot labels 452 if labels: 453 if 100 > len(lons) > 1: 454 for name, xpt, ypt, _colorpt in zip(labels, x, y, color): 455 # Check if the point can actually be seen with the current bmap 456 # projection. The bmap object will set the coordinates to very 457 # large values if it cannot project a point. 458 if xpt > 1e25: 459 continue 460 ax.text(xpt, ypt, name, weight="heavy", 461 color="k", zorder=100, 462 path_effects=[ 463 patheffects.withStroke(linewidth=3, 464 foreground="white")]) 465 elif len(lons) == 1: 466 ax.text(x[0], y[0], labels[0], weight="heavy", color="k", 467 path_effects=[ 468 patheffects.withStroke(linewidth=3, 469 foreground="white")]) 470 471 # scatter plot is removing valid x/y points with invalid color value, 472 # so we plot those points separately. 473 try: 474 nan_points = np.isnan(np.array(color, dtype=np.float)) 475 except ValueError: 476 # `color' was not a list of values, but a list of colors. 477 pass 478 else: 479 if nan_points.any(): 480 x_ = np.array(x)[nan_points] 481 y_ = np.array(y)[nan_points] 482 size_ = np.array(size)[nan_points] 483 bmap.scatter(x_, y_, marker=marker, s=size_, c="0.3", 484 zorder=10, cmap=None) 485 scatter = bmap.scatter(x, y, marker=marker, s=size, c=color, zorder=10, 486 cmap=colormap) 487 488 if title: 489 ax.set_title(title) 490 491 return scatter 492 493 494def plot_cartopy(lons, lats, size, color, labels=None, projection='global', 495 resolution='110m', continent_fill_color='0.8', 496 water_fill_color='1.0', colormap=None, colorbar=None, 497 marker="o", title=None, colorbar_ticklabel_format=None, 498 show=True, proj_kwargs=None, **kwargs): # @UnusedVariable 499 """ 500 Creates a Cartopy plot with a data point scatter plot. 501 502 :type lons: list/tuple of floats 503 :param lons: Longitudes of the data points. 504 :type lats: list/tuple of floats 505 :param lats: Latitudes of the data points. 506 :type size: float or list/tuple of floats 507 :param size: Size of the individual points in the scatter plot. 508 :type color: list/tuple of floats (or objects that can be 509 converted to floats, like e.g. 510 :class:`~obspy.core.utcdatetime.UTCDateTime`) 511 :param color: Color information of the individual data points to be 512 used in the specified color map (e.g. origin depths, 513 origin times). 514 :type labels: list/tuple of str 515 :param labels: Annotations for the individual data points. 516 :type projection: str, optional 517 :param projection: The map projection. 518 Currently supported are: 519 520 * ``"global"`` (Will plot the whole world using 521 :class:`~cartopy.crs.Mollweide`.) 522 * ``"ortho"`` (Will center around the mean lat/long using 523 :class:`~cartopy.crs.Orthographic`.) 524 * ``"local"`` (Will plot around local events using 525 :class:`~cartopy.crs.AlbersEqualArea`.) 526 * Any other Cartopy :class:`~cartopy.crs.Projection`. An instance 527 of this class will be created using the supplied ``proj_kwargs``. 528 529 Defaults to "global" 530 :type resolution: str, optional 531 :param resolution: Resolution of the boundary database to use. Will be 532 passed directly to the Cartopy module. Possible values are: 533 534 * ``"110m"`` 535 * ``"50m"`` 536 * ``"10m"`` 537 538 Defaults to ``"110m"``. For compatibility, you may also specify any of 539 the Basemap resolutions defined in :func:`plot_basemap`. 540 :type continent_fill_color: Valid matplotlib color, optional 541 :param continent_fill_color: Color of the continents. Defaults to 542 ``"0.9"`` which is a light gray. 543 :type water_fill_color: Valid matplotlib color, optional 544 :param water_fill_color: Color of all water bodies. 545 Defaults to ``"white"``. 546 :type colormap: str, any matplotlib colormap, optional 547 :param colormap: The colormap for color-coding the events as provided 548 in `color` kwarg. 549 The event with the smallest `color` property will have the 550 color of one end of the colormap and the event with the highest 551 `color` property the color of the other end with all other events 552 in between. 553 Defaults to None which will use the default matplotlib colormap. 554 :type colorbar: bool, optional 555 :param colorbar: When left `None`, a colorbar is plotted if more than one 556 object is plotted. Using `True`/`False` the colorbar can be forced 557 on/off. 558 :type title: str 559 :param title: Title above plot. 560 :type colorbar_ticklabel_format: str or function or 561 subclass of :class:`matplotlib.ticker.Formatter` 562 :param colorbar_ticklabel_format: Format string or Formatter used to format 563 colorbar tick labels. 564 :type show: bool 565 :param show: Whether to show the figure after plotting or not. Can be used 566 to do further customization of the plot before showing it. 567 :type proj_kwargs: dict 568 :param proj_kwargs: Keyword arguments to pass to the Cartopy 569 :class:`~cartopy.ccrs.Projection`. In this dictionary, you may specify 570 ``central_longitude='auto'`` or ``central_latitude='auto'`` to have 571 this function calculate the latitude or longitude as it would for other 572 projections. Some arguments may be ignored if you choose one of the 573 built-in ``projection`` choices. 574 """ 575 import matplotlib.pyplot as plt 576 577 if isinstance(color[0], (datetime.datetime, UTCDateTime)): 578 datetimeplot = True 579 color = [date2num(getattr(t, 'datetime', t)) for t in color] 580 else: 581 datetimeplot = False 582 583 fig = plt.figure() 584 585 # The colorbar should only be plotted if more then one event is 586 # present. 587 if colorbar is not None: 588 show_colorbar = colorbar 589 else: 590 if len(lons) > 1 and hasattr(color, "__len__") and \ 591 not isinstance(color, (str, native_str)): 592 show_colorbar = True 593 else: 594 show_colorbar = False 595 596 if projection == "local": 597 ax_x0, ax_width = 0.10, 0.80 598 elif projection == "global": 599 ax_x0, ax_width = 0.01, 0.98 600 else: 601 ax_x0, ax_width = 0.05, 0.90 602 603 proj_kwargs = proj_kwargs or {} 604 if projection == 'global': 605 proj_kwargs['central_longitude'] = np.mean(lons) 606 proj = ccrs.Mollweide(**proj_kwargs) 607 elif projection == 'ortho': 608 proj_kwargs['central_latitude'] = np.mean(lats) 609 proj_kwargs['central_longitude'] = mean_longitude(lons) 610 proj = ccrs.Orthographic(**proj_kwargs) 611 elif projection == 'local': 612 if min(lons) < -150 and max(lons) > 150: 613 max_lons = max(np.array(lons) % 360) 614 min_lons = min(np.array(lons) % 360) 615 else: 616 max_lons = max(lons) 617 min_lons = min(lons) 618 lat_0 = max(lats) / 2. + min(lats) / 2. 619 lon_0 = max_lons / 2. + min_lons / 2. 620 if lon_0 > 180: 621 lon_0 -= 360 622 deg2m_lat = 2 * np.pi * 6371 * 1000 / 360 623 deg2m_lon = deg2m_lat * np.cos(lat_0 / 180 * np.pi) 624 if len(lats) > 1: 625 height = (max(lats) - min(lats)) * deg2m_lat 626 width = (max_lons - min_lons) * deg2m_lon 627 margin = 0.2 * (width + height) 628 height += margin 629 width += margin 630 else: 631 height = 2.0 * deg2m_lat 632 width = 5.0 * deg2m_lon 633 # Do intelligent aspect calculation for local projection 634 # adjust to figure dimensions 635 w, h = fig.get_size_inches() 636 aspect = w / h 637 if show_colorbar: 638 aspect *= 1.2 639 if width / height < aspect: 640 width = height * aspect 641 else: 642 height = width / aspect 643 644 proj_kwargs['central_latitude'] = lat_0 645 proj_kwargs['central_longitude'] = lon_0 646 proj_kwargs['standard_parallels'] = [lat_0, lat_0] 647 proj = ccrs.AlbersEqualArea(**proj_kwargs) 648 649 # User-supplied projection. 650 elif isinstance(projection, type): 651 if 'central_longitude' in proj_kwargs: 652 if proj_kwargs['central_longitude'] == 'auto': 653 proj_kwargs['central_longitude'] = mean_longitude(lons) 654 if 'central_latitude' in proj_kwargs: 655 if proj_kwargs['central_latitude'] == 'auto': 656 proj_kwargs['central_latitude'] = np.mean(lats) 657 if 'pole_longitude' in proj_kwargs: 658 if proj_kwargs['pole_longitude'] == 'auto': 659 proj_kwargs['pole_longitude'] = np.mean(lons) 660 if 'pole_latitude' in proj_kwargs: 661 if proj_kwargs['pole_latitude'] == 'auto': 662 proj_kwargs['pole_latitude'] = np.mean(lats) 663 664 proj = projection(**proj_kwargs) 665 666 else: 667 msg = "Projection '%s' not supported." % projection 668 raise ValueError(msg) 669 670 if show_colorbar: 671 map_ax = fig.add_axes([ax_x0, 0.13, ax_width, 0.77], projection=proj) 672 cm_ax = fig.add_axes([ax_x0, 0.05, ax_width, 0.05]) 673 plt.sca(map_ax) 674 else: 675 ax_y0, ax_height = 0.05, 0.85 676 if projection == "local": 677 ax_y0 += 0.05 678 ax_height -= 0.05 679 map_ax = fig.add_axes([ax_x0, ax_y0, ax_width, ax_height], 680 projection=proj) 681 682 if projection == 'local': 683 x0, y0 = proj.transform_point(lon_0, lat_0, proj.as_geodetic()) 684 map_ax.set_xlim(x0 - width / 2, x0 + width / 2) 685 map_ax.set_ylim(y0 - height / 2, y0 + height / 2) 686 else: 687 map_ax.set_global() 688 689 # Pick features at specified resolution. 690 resolution = _CARTOPY_RESOLUTIONS[resolution] 691 try: 692 borders, land, ocean = _CARTOPY_FEATURES[resolution] 693 except KeyError: 694 borders = cfeature.NaturalEarthFeature(cfeature.BORDERS.category, 695 cfeature.BORDERS.name, 696 resolution, 697 edgecolor='none', 698 facecolor='none') 699 land = cfeature.NaturalEarthFeature(cfeature.LAND.category, 700 cfeature.LAND.name, resolution, 701 edgecolor='face', facecolor='none') 702 ocean = cfeature.NaturalEarthFeature(cfeature.OCEAN.category, 703 cfeature.OCEAN.name, resolution, 704 edgecolor='face', 705 facecolor='none') 706 _CARTOPY_FEATURES[resolution] = (borders, land, ocean) 707 708 # Draw coast lines, country boundaries, fill continents. 709 if MATPLOTLIB_VERSION >= [2, 0, 0]: 710 map_ax.set_facecolor(water_fill_color) 711 else: 712 map_ax.set_axis_bgcolor(water_fill_color) 713 map_ax.add_feature(ocean, facecolor=water_fill_color) 714 map_ax.add_feature(land, facecolor=continent_fill_color) 715 map_ax.add_feature(borders, edgecolor='0.75') 716 map_ax.coastlines(resolution=resolution, color='0.4') 717 718 # Draw grid lines - TODO: draw_labels=True doesn't work yet. 719 if projection == 'local': 720 map_ax.gridlines() 721 else: 722 # Draw lat/lon grid lines every 30 degrees. 723 map_ax.gridlines(xlocs=range(-180, 181, 30), ylocs=range(-90, 91, 30)) 724 725 # Plot labels 726 if labels and len(lons) > 0: 727 with map_ax.hold_limits(): 728 for name, xpt, ypt, _colorpt in zip(labels, lons, lats, color): 729 map_ax.text(xpt, ypt, name, weight="heavy", color="k", 730 zorder=100, transform=ccrs.PlateCarree(), 731 path_effects=[ 732 patheffects.withStroke(linewidth=3, 733 foreground="white")]) 734 735 scatter = map_ax.scatter(lons, lats, marker=marker, s=size, c=color, 736 zorder=10, cmap=colormap, 737 transform=ccrs.PlateCarree()) 738 739 if title: 740 plt.suptitle(title) 741 742 # Only show the colorbar for more than one event. 743 if show_colorbar: 744 if colorbar_ticklabel_format is not None: 745 if isinstance(colorbar_ticklabel_format, (str, native_str)): 746 formatter = FormatStrFormatter(colorbar_ticklabel_format) 747 elif hasattr(colorbar_ticklabel_format, '__call__'): 748 formatter = FuncFormatter(colorbar_ticklabel_format) 749 elif isinstance(colorbar_ticklabel_format, Formatter): 750 formatter = colorbar_ticklabel_format 751 locator = MaxNLocator(5) 752 else: 753 if datetimeplot: 754 locator = AutoDateLocator() 755 formatter = AutoDateFormatter(locator) 756 # Compat with old matplotlib versions. 757 if hasattr(formatter, "scaled"): 758 formatter.scaled[1 / (24. * 60.)] = '%H:%M:%S' 759 else: 760 locator = None 761 formatter = None 762 cb = Colorbar(cm_ax, scatter, cmap=colormap, 763 orientation='horizontal', 764 ticks=locator, 765 format=formatter) 766 # Compat with old matplotlib versions. 767 if hasattr(cb, "update_ticks"): 768 cb.update_ticks() 769 770 if show: 771 plt.show() 772 773 return fig 774 775 776def plot_map(method, *args, **kwargs): 777 ''' 778 Creates a map plot with a data point scatter plot. 779 780 :type method: str 781 :param method: Method to use for plotting. Possible values are: 782 783 * ``'basemap'`` to use the Basemap library. For other arguments, see 784 the :func:`plot_basemap` function. 785 * ``'cartopy'`` to use the Cartopy library. For other arguments, see 786 the :func:`plot_cartopy` function. 787 * ``None`` to use either the Basemap or Cartopy library, whichever is 788 available. 789 ''' 790 791 if method is None: 792 if HAS_BASEMAP: 793 return plot_basemap(*args, **kwargs) 794 elif HAS_CARTOPY: 795 return plot_cartopy(*args, **kwargs) 796 else: 797 raise ImportError('Neither Basemap nor Cartopy could be imported.') 798 elif method == 'basemap': 799 if not HAS_BASEMAP: 800 raise ImportError('Basemap cannot be imported but was explicitly ' 801 'requested.') 802 return plot_basemap(*args, **kwargs) 803 elif method == 'cartopy': 804 if not HAS_CARTOPY: 805 raise ImportError('Cartopy cannot be imported but was explicitly ' 806 'requested.') 807 return plot_cartopy(*args, **kwargs) 808 else: 809 raise ValueError("The method argument must be either 'basemap' or " 810 "'cartopy', not '%s'." % (method, )) 811