From e5e4a2e26e83f78d9998e4321882b2f8a22c17f7 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sat, 21 Nov 2020 11:27:57 +0100 Subject: [PATCH] Deprecate Basemap dependency (#180) * Deprecate Basemap dependency * Plot coastline * Support AzimuthalEquidistant projection * Catch exception when Cartopy is not installed * Add option to draw lat lon labels * Support Geostationary projection * Fix tests * Fix SyntaxWarning * Refactor kwargs options * Add plot_map optional argument * Fix interface to basemaps (Cartopy is now used by default) --- examples/my_first_nowcast.ipynb | 6 +- pysteps/nowcasts/sseps.py | 2 +- pysteps/tests/test_plt_cartopy.py | 37 ++-- pysteps/tests/test_plt_motionfields.py | 70 ++---- pysteps/tests/test_plt_precipfields.py | 18 +- pysteps/visualization/animations.py | 108 ++++++---- pysteps/visualization/basemaps.py | 283 +++++++++---------------- pysteps/visualization/motionfields.py | 60 +++--- pysteps/visualization/precipfields.py | 65 +++--- pysteps/visualization/utils.py | 64 +----- 10 files changed, 285 insertions(+), 428 deletions(-) diff --git a/examples/my_first_nowcast.ipynb b/examples/my_first_nowcast.ipynb index 338ad1f50..3113c4977 100644 --- a/examples/my_first_nowcast.ipynb +++ b/examples/my_first_nowcast.ipynb @@ -417,7 +417,7 @@ "\n", "# Plot the last rainfall field in the \"training\" data.\n", "# train_precip[-1] -> Last available composite for nowcasting.\n", - "plot_precip_field(train_precip[-1], geodata=metadata, axis=\"off\", map=\"cartopy\")\n", + "plot_precip_field(train_precip[-1], geodata=metadata, axis=\"off\")\n", "plt.show() # (This line is actually not needed if you are using jupyter notebooks)" ] }, @@ -699,7 +699,7 @@ "source": [ "# Plot precipitation at the end of the forecast period.\n", "plt.figure(figsize=(9, 5), dpi=100)\n", - "plot_precip_field(precip_forecast[-1], geodata=metadata, axis=\"off\", map='cartopy')\n", + "plot_precip_field(precip_forecast[-1], geodata=metadata, axis=\"off\")\n", "plt.show()" ] }, @@ -829,4 +829,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/pysteps/nowcasts/sseps.py b/pysteps/nowcasts/sseps.py index 1db76fe2b..e3c03ed4f 100644 --- a/pysteps/nowcasts/sseps.py +++ b/pysteps/nowcasts/sseps.py @@ -271,7 +271,7 @@ def forecast( print("dask imported: %s" % ("yes" if dask_imported else "no")) print("num workers: %d" % num_workers) - if vel_pert_method is "bps": + if vel_pert_method == "bps": vp_par = vel_pert_kwargs.get( "p_pert_par", noise.motion.get_default_params_bps_par() ) diff --git a/pysteps/tests/test_plt_cartopy.py b/pysteps/tests/test_plt_cartopy.py index 7ad0ce2ce..ea3eacd29 100644 --- a/pysteps/tests/test_plt_cartopy.py +++ b/pysteps/tests/test_plt_cartopy.py @@ -5,43 +5,42 @@ from pysteps.visualization import plot_precip_field from pysteps.utils import to_rainrate from pysteps.tests.helpers import get_precipitation_fields -import matplotlib.pyplot as pl +import matplotlib.pyplot as plt -pytest.importorskip("cartopy") -plt_arg_names = ("source", "plot_map", "drawlonlatlines", "lw") +plt_arg_names = ("source", "map_kwargs", "pass_geodata") plt_arg_values = [ - ("mch", "cartopy", False, 0.5), - ("mch", "cartopy", True, 1.0), - ("bom", "cartopy", True, 0.5), - ("fmi", "cartopy", True, 0.5), - ("knmi", "cartopy", True, 0.5), - ("opera", "cartopy", True, 0.5), - ("mrms", "cartopy", True, 0.5), + ("mch", {"drawlonlatlines": False, "lw": 0.5, "plot_map": None}, False), + ("mch", {"drawlonlatlines": False, "lw": 0.5, "plot_map": "cartopy"}, False), + ("mch", {"drawlonlatlines": False, "lw": 0.5}, True), + ("mch", {"drawlonlatlines": True, "lw": 1.0}, True), + ("bom", {"drawlonlatlines": True, "lw": 0.5}, True), + ("fmi", {"drawlonlatlines": True, "lw": 0.5}, True), + ("knmi", {"drawlonlatlines": True, "lw": 0.5}, True), + ("opera", {"drawlonlatlines": True, "lw": 0.5}, True), + ("mrms", {"drawlonlatlines": True, "lw": 0.5}, True), + ("saf", {"drawlonlatlines": True, "lw": 0.5}, True), ] @pytest.mark.parametrize(plt_arg_names, plt_arg_values) -def test_visualization_plot_precip_field(source, plot_map, drawlonlatlines, lw): +def test_visualization_plot_precip_field(source, map_kwargs, pass_geodata): field, metadata = get_precipitation_fields(0, 0, True, True, None, source) field = field.squeeze() field, __ = to_rainrate(field, metadata) + if not pass_geodata: + metadata = None + ax = plot_precip_field( - field, - type="intensity", - geodata=metadata, - plot_map=plot_map, - drawlonlatlines=drawlonlatlines, - lw=lw, + field, type="intensity", geodata=metadata, map_kwargs=map_kwargs, ) - pl.close() if __name__ == "__main__": for i, args in enumerate(plt_arg_values): test_visualization_plot_precip_field(*args) - pl.show() + plt.show() diff --git a/pysteps/tests/test_plt_motionfields.py b/pysteps/tests/test_plt_motionfields.py index a94e3cbb9..b4e06dfe3 100644 --- a/pysteps/tests/test_plt_motionfields.py +++ b/pysteps/tests/test_plt_motionfields.py @@ -14,34 +14,29 @@ "axis", "step", "quiver_kwargs", - "plot_map", - "drawlonlatlines", - "lw", + "map_kwargs", "upscale", + "pass_geodata", ) arg_values_quiver = [ - (None, "off", 10, None, None, False, 0.5, None), - ("bom", "on", 10, None, None, False, 0.5, 4000), - ("bom", "on", 10, None, "cartopy", True, 0.5, 4000), - ("mch", "on", 20, None, "cartopy", False, 0.5, 2000), - ("bom", "off", 10, None, "basemap", False, 0.5, 4000), + (None, "off", 10, {}, {"drawlonlatlines": False, "lw": 0.5}, None, False), + ("bom", "on", 10, {}, {"drawlonlatlines": False, "lw": 0.5}, 4000, False), + ("bom", "on", 10, {}, {"drawlonlatlines": True, "lw": 0.5}, 4000, True), + ("mch", "on", 20, {}, {"drawlonlatlines": False, "lw": 0.5}, 2000, True), ] @pytest.mark.parametrize(arg_names_quiver, arg_values_quiver) def test_visualization_motionfields_quiver( - source, axis, step, quiver_kwargs, plot_map, drawlonlatlines, lw, upscale, + source, axis, step, quiver_kwargs, map_kwargs, upscale, pass_geodata, ): - if plot_map == "cartopy": - pytest.importorskip("cartopy") - elif plot_map == "basemap": - pytest.importorskip("basemap") - if source is not None: fields, geodata = get_precipitation_fields(0, 2, False, True, upscale, source) - ax = plot_precip_field(fields[-1], geodata=geodata, plot_map=plot_map,) + if not pass_geodata: + geodata = None + ax = plot_precip_field(fields[-1], geodata=geodata) oflow_method = motion.get_method("LK") UV = oflow_method(fields) @@ -54,49 +49,35 @@ def test_visualization_motionfields_quiver( U, V = np.meshgrid(u, v) UV = np.concatenate([U[None, :], V[None, :]]) - __ = quiver( - UV, - ax, - geodata, - axis, - step, - quiver_kwargs, - plot_map=plot_map, - drawlonlatlines=drawlonlatlines, - lw=lw, - ) + __ = quiver(UV, ax, geodata, axis, step, quiver_kwargs, map_kwargs=map_kwargs,) arg_names_streamplot = ( "source", "axis", "streamplot_kwargs", - "plot_map", - "drawlonlatlines", - "lw", + "map_kwargs", "upscale", + "pass_geodata", ) arg_values_streamplot = [ - (None, "off", None, None, False, 0.5, None), - ("bom", "on", None, None, False, 0.5, 4000), - ("bom", "on", {"density": 0.1}, "cartopy", True, 0.5, 4000), + (None, "off", {}, {"drawlonlatlines": False, "lw": 0.5}, None, False), + ("bom", "on", {}, {"drawlonlatlines": False, "lw": 0.5}, 4000, False), + ("bom", "on", {"density": 0.5}, {"drawlonlatlines": True, "lw": 0.5}, 4000, True), ] @pytest.mark.parametrize(arg_names_streamplot, arg_values_streamplot) def test_visualization_motionfields_streamplot( - source, axis, streamplot_kwargs, plot_map, drawlonlatlines, lw, upscale, + source, axis, streamplot_kwargs, map_kwargs, upscale, pass_geodata ): - if plot_map == "cartopy": - pytest.importorskip("cartopy") - elif plot_map == "basemap": - pytest.importorskip("basemap") - if source is not None: fields, geodata = get_precipitation_fields(0, 2, False, True, upscale, source) - ax = plot_precip_field(fields[-1], geodata=geodata, plot_map=plot_map,) + if not pass_geodata: + pass_geodata = None + ax = plot_precip_field(fields[-1], geodata=geodata) oflow_method = motion.get_method("LK") UV = oflow_method(fields) @@ -109,16 +90,7 @@ def test_visualization_motionfields_streamplot( U, V = np.meshgrid(u, v) UV = np.concatenate([U[None, :], V[None, :]]) - __ = streamplot( - UV, - ax, - geodata, - axis, - streamplot_kwargs, - plot_map=plot_map, - drawlonlatlines=drawlonlatlines, - lw=lw, - ) + __ = streamplot(UV, ax, geodata, axis, streamplot_kwargs, map_kwargs=map_kwargs,) if __name__ == "__main__": diff --git a/pysteps/tests/test_plt_precipfields.py b/pysteps/tests/test_plt_precipfields.py index 9c58756d8..7d49bf1a7 100644 --- a/pysteps/tests/test_plt_precipfields.py +++ b/pysteps/tests/test_plt_precipfields.py @@ -6,7 +6,7 @@ from pysteps.utils import conversion from pysteps.postprocessing import ensemblestats from pysteps.tests.helpers import get_precipitation_fields -import matplotlib.pyplot as pl +import matplotlib.pyplot as plt plt_arg_names = ( "source", @@ -28,16 +28,7 @@ ("bom", "intensity", None, "pysteps", None, None, True, "on"), ("fmi", "intensity", None, "pysteps", None, None, True, "on"), ("knmi", "intensity", None, "pysteps", None, None, True, "on"), - ( - "knmi", - "intensity", - [2e2, -4.1e3, 5e2, -3.8e3], - "pysteps", - None, - None, - True, - "on", - ), + ("knmi", "intensity", [300, 300, 500, 500], "pysteps", None, None, True, "on",), ("opera", "intensity", None, "pysteps", None, None, True, "on"), ("saf", "intensity", None, "pysteps", None, None, True, "on"), ] @@ -70,7 +61,7 @@ def test_visualization_plot_precip_field( field, type=type, bbox=bbox, - geodata=metadata, + geodata=None, colorscale=colorscale, probthr=probthr, units=metadata["unit"], @@ -78,11 +69,10 @@ def test_visualization_plot_precip_field( colorbar=colorbar, axis=axis, ) - pl.close() if __name__ == "__main__": for i, args in enumerate(plt_arg_values): test_visualization_plot_precip_field(*args) - pl.show() + plt.show() diff --git a/pysteps/visualization/animations.py b/pysteps/visualization/animations.py index e02679f13..a831f409d 100644 --- a/pysteps/visualization/animations.py +++ b/pysteps/visualization/animations.py @@ -10,6 +10,8 @@ animate """ +import os + import matplotlib.pylab as plt import numpy as np import pysteps as st @@ -35,7 +37,8 @@ def animate( fig_dpi=150, fig_format="png", path_outputs="", - **kwargs, + motion_kwargs={}, + map_kwargs={}, ): """Function to animate observations and forecasts in pysteps. @@ -58,7 +61,7 @@ def animate( Optional, the motion field used for the forecast. motion_plot : string The method to plot the motion field. - geodata : dictionary + geodata : dictionary or None Optional dictionary containing geographical information about the field. If geodata is not None, it must contain the following key-value pairs: @@ -113,10 +116,12 @@ def animate( matplotlib.pyplot.savefig. Applicable if savefig is True. path_outputs : string Path to folder where to save the frames. - kwargs : dict + motion_kwargs : dict + Optional keyword arguments that are supplied to + :py:func:`pysteps.visualization.precipfields.plot_precip_field` or + :py:func:`pysteps.visualization.motionfields.quiver`. + map_kwargs : dict Optional keyword arguments that are supplied to - :py:func:`pysteps.visualization.precipfields.plot_precip_field`, - :py:func:`pysteps.visualization.motionfields.quiver`, and :py:func:`pysteps.visualization.motionfields.streamplot`. Returns @@ -176,7 +181,7 @@ def animate( units=units, probthr=prob_thr, title=title, - **kwargs, + map_kwargs=map_kwargs, ) else: title += "Observed Rainfall" @@ -187,30 +192,38 @@ def animate( colorscale=colorscale, title=title, colorbar=colorbar, - **kwargs, + map_kwargs=map_kwargs, ) if UV is not None and motion_plot is not None: if motion_plot.lower() == "quiver": - st.plt.quiver(UV, ax=ax, geodata=geodata, **kwargs) + st.plt.quiver(UV, ax=ax, geodata=geodata, **motion_kwargs) elif motion_plot.lower() == "streamplot": - st.plt.streamplot(UV, ax=ax, geodata=geodata, **kwargs) + st.plt.streamplot( + UV, ax=ax, geodata=geodata, **motion_kwargs + ) if savefig & (loop == 0): if type == "prob": - figname = "%s/%s_frame_%02d_binmap_%.1f.%s" % ( - path_outputs, - startdate_str, - i, - prob_thr, - fig_format, + figname = os.path.join( + "%s, %s_frame_%02d_binmap_%.1f.%s" + % ( + path_outputs, + startdate_str, + i, + prob_thr, + fig_format, + ) ) else: - figname = "%s/%s_frame_%02d.%s" % ( - path_outputs, - startdate_str, - i, - fig_format, + figname = os.path.join( + "%s, %s_frame_%02d.%s" + % ( + path_outputs, + startdate_str, + i, + fig_format, + ) ) plt.savefig(figname, bbox_inches="tight", dpi=fig_dpi) print(figname, "saved.") @@ -236,7 +249,7 @@ def animate( units=units, probthr=prob_thr, title=title, - **kwargs, + map_kwargs=map_kwargs, ) elif type == "mean": title += "Forecast Ensemble Mean" @@ -250,7 +263,7 @@ def animate( title=title, colorscale=colorscale, colorbar=colorbar, - **kwargs, + map_kwargs=map_kwargs, ) else: title += "Forecast Rainfall" @@ -261,14 +274,16 @@ def animate( title=title, colorscale=colorscale, colorbar=colorbar, - **kwargs, + map_kwargs=map_kwargs, ) if UV is not None and motion_plot is not None: if motion_plot.lower() == "quiver": - st.plt.quiver(UV, ax=ax, geodata=geodata, **kwargs) + st.plt.quiver(UV, ax=ax, geodata=geodata, **motion_kwargs) elif motion_plot.lower() == "streamplot": - st.plt.streamplot(UV, ax=ax, geodata=geodata, **kwargs) + st.plt.streamplot( + UV, ax=ax, geodata=geodata, **motion_kwargs + ) if leadtime is not None: plt.text( @@ -291,27 +306,36 @@ def animate( if savefig & (loop == 0): if type == "prob": - figname = "%s/%s_frame_%02d_probmap_%.1f.%s" % ( - path_outputs, - startdate_str, - i, - prob_thr, - fig_format, + figname = os.path.join( + "%s, %s_frame_%02d_probmap_%.1f.%s" + % ( + path_outputs, + startdate_str, + i, + prob_thr, + fig_format, + ) ) elif type == "mean": - figname = "%s/%s_frame_%02d_ensmean.%s" % ( - path_outputs, - startdate_str, - i, - fig_format, + figname = os.path.join( + "%s, %s_frame_%02d_ensmean.%s" + % ( + path_outputs, + startdate_str, + i, + fig_format, + ) ) else: - figname = "%s/%s_member_%02d_frame_%02d.%s" % ( - path_outputs, - startdate_str, - (n + 1), - i, - fig_format, + figname = os.path.join( + "%s, %s_member_%02d_frame_%02d.%s" + % ( + path_outputs, + startdate_str, + (n + 1), + i, + fig_format, + ) ) plt.savefig(figname, bbox_inches="tight", dpi=fig_dpi) print(figname, "saved.") diff --git a/pysteps/visualization/basemaps.py b/pysteps/visualization/basemaps.py index f3cb3eb45..507b63f7a 100644 --- a/pysteps/visualization/basemaps.py +++ b/pysteps/visualization/basemaps.py @@ -2,13 +2,12 @@ pysteps.visualization.basemaps ============================== -Methods for plotting geographical maps using Cartopy or Basemap. +Methods for plotting geographical maps using Cartopy. .. autosummary:: :toctree: ../generated/ plot_geography - plot_map_basemap plot_map_cartopy """ from matplotlib import gridspec @@ -17,253 +16,151 @@ import warnings from pysteps.exceptions import MissingOptionalDependency -try: - from mpl_toolkits.basemap import Basemap - basemap_imported = True -except ImportError: - basemap_imported = False try: import cartopy.crs as ccrs import cartopy.feature as cfeature - cartopy_imported = True + CARTOPY_IMPORTED = True except ImportError: - cartopy_imported = False + CARTOPY_IMPORTED = False try: import pyproj - pyproj_imported = True + PYPROJ_IMPORTED = True except ImportError: - pyproj_imported = False + PYPROJ_IMPORTED = False from . import utils +VALID_BASEMAPS = ["cartopy"] + + def plot_geography( - proj4str, extent, plot_map="cartopy", lw=0.5, drawlonlatlines=False, **kwargs + proj4str, extent, lw=0.5, drawlonlatlines=False, drawlonlatlabels=True, **kwargs ): """ - Plot geographical map using either cartopy_ or basemap_ in a chosen projection. + Plot geographical map using cartopy_ in a chosen projection. .. _cartopy: https://scitools.org.uk/cartopy/docs/latest - .. _basemap: https://matplotlib.org/basemap - .. _SubplotSpec: https://matplotlib.org/api/_as_gen/matplotlib.gridspec.SubplotSpec.html Parameters ---------- - proj4str : str + proj4str: str The PROJ.4-compatible projection string. extent: scalars (left, right, bottom, top) The bounding box in proj4str coordinates. - plot_map : {'cartopy', 'basemap'}, optional - The type of basemap, either 'cartopy_' or 'basemap_'. lw: float, optional Linewidth of the map (administrative boundaries and coastlines). - drawlonlatlines : bool, optional - If set to True, draw longitude and latitude lines. Applicable if plot_map is - 'basemap' or 'cartopy'. + drawlonlatlines: bool, optional + If set to True, draw longitude and latitude lines. + drawlonlatlabels: bool, optional + If set to True, draw longitude and latitude labels. Valid only if + 'drawlonlatlines' is True. Other parameters ---------------- - resolution : str, optional - The resolution of the map, see the documentation of `basemap`_. - Applicable if 'plot_map' is 'basemap'. Default ``'l'``. - scale_args : list, optional - If not None, a map scale bar is drawn with basemap_scale_args supplied - to mpl_toolkits.basemap.Basemap.drawmapscale. Applicable if 'plot_map' is - 'basemap'. Default ``None``. - scale : {'10m', '50m', '110m'}, optional - The scale (resolution) of the plot_map. The available options are '10m', - '50m', and '110m'. Applicable if plot_map is 'cartopy'. Default ``'50m'`` - subplot : tuple or SubplotSpec_ instance, optional - The cartopy subplot to plot into. Applicable if plot_map is 'cartopy'. - Default ``'(1, 1, 1)'`` + plot_map : {'cartopy', None}, optional + The type of basemap, either 'cartopy_' or None. If None, the figure + axis is returned without any basemap drawn. Default ``'cartopy'``. + scale: {'10m', '50m', '110m'}, optional + The scale (resolution). Applicable if 'plot_map' is 'cartopy'. + The available options are '10m', '50m', and '110m'. Default ``'50m'``. + subplot: tuple or SubplotSpec_ instance, optional + The cartopy subplot to plot into. Applicable if 'plot_map' is 'cartopy'. + Default ``'(1, 1, 1)'``. Returns ------- - ax : fig Axes_ - Cartopy or Basemap axes. - regular_grid : bool - Whether the projection allows plotting a regular grid. - Returns False in case a fall-back projection is used. + ax: fig Axes_ + Cartopy axes. """ - if plot_map not in ["basemap", "cartopy"]: + plot_map = kwargs.get("plot_map", "cartopy") + + if plot_map is None: + return plt.axes() + + if plot_map not in VALID_BASEMAPS: raise ValueError( - "unknown plot_map method %s: must be" + " 'basemap' or 'cartopy'" % plot_map - ) - if plot_map == "basemap" and not basemap_imported: - raise MissingOptionalDependency( - "plot_map='basemap' option passed to plot_geography function " - "but the basemap package is not installed" - ) - if plot_map == "cartopy" and not cartopy_imported: - raise MissingOptionalDependency( - "plot_map='cartopy' option passed to plot_geography function " - "but the cartopy package is not installed" - ) - if plot_map is not None and not pyproj_imported: - raise MissingOptionalDependency( - "plot_map!=None option passed to plot_geography function " - "but the pyproj package is not installed" + f"unsupported plot_map method {plot_map}. Supported basemaps: " + f"{VALID_BASEMAPS}" ) - if plot_map == "basemap": - basemap_resolution = kwargs.get("resolution", "l") - basemap_scale_args = kwargs.get("scale_args", None) - pr = pyproj.Proj(proj4str) - x1, x2, y1, y2 = extent[0], extent[1], extent[2], extent[3] - ll_lon, ll_lat = pr(x1, y1, inverse=True) - ur_lon, ur_lat = pr(x2, y2, inverse=True) - - bm_params = utils.proj4_to_basemap(proj4str) - bm_params["llcrnrlon"] = ll_lon - bm_params["llcrnrlat"] = ll_lat - bm_params["urcrnrlon"] = ur_lon - bm_params["urcrnrlat"] = ur_lat - bm_params["resolution"] = basemap_resolution - - ax = plot_map_basemap(bm_params, drawlonlatlines=drawlonlatlines, lw=lw) - - if basemap_scale_args is not None: - ax.drawmapscale(*basemap_scale_args, fontsize=6, yoffset=10000) - else: - cartopy_scale = kwargs.get("scale", "50m") - cartopy_subplot = kwargs.get("subplot", (1, 1, 1)) - crs = utils.proj4_to_cartopy(proj4str) - - ax = plot_map_cartopy( - crs, - extent, - cartopy_scale, - drawlonlatlines=drawlonlatlines, - lw=lw, - subplot=cartopy_subplot, + if plot_map == "cartopy" and not CARTOPY_IMPORTED: + raise MissingOptionalDependency( + "the cartopy package is required to plot the geographical map " + " but it is not installed" ) - return ax - - -def plot_map_basemap( - bm_params, - drawlonlatlines=False, - coastlinecolor=(0.3, 0.3, 0.3), - countrycolor=(0.3, 0.3, 0.3), - continentcolor=(0.95, 0.95, 0.85), - lakecolor=(0.65, 0.75, 0.9), - rivercolor=(0.65, 0.75, 0.9), - mapboundarycolor=(0.65, 0.75, 0.9), - lw=0.5, -): - """ - Plot coastlines, countries, rivers and meridians/parallels using Basemap. - - Parameters - ---------- - bm_params : optional - Optional arguments for the Basemap class constructor: - https://basemaptutorial.readthedocs.io/en/latest/basemap.html - drawlonlatlines : bool - Whether to plot longitudes and latitudes. - coastlinecolor : scalars (r, g, b) - Coastline color. - countrycolor : scalars (r, g, b) - Countrycolor color. - continentcolor : scalars (r, g, b) - Continentcolor color. - lakecolor : scalars (r, g, b) - Lakecolor color. - rivercolor : scalars (r, g, b) - Rivercolor color. - mapboundarycolor : scalars (r, g, b) - Mapboundarycolor color. - lw : float - Line width. - - Returns - ------- - ax : axes - Basemap axes. - """ - if not basemap_imported: + if not PYPROJ_IMPORTED: raise MissingOptionalDependency( - "plot_map='basemap' option passed to plot_map_basemap function " - "but the basemap package is not installed" + "the pyproj package is required to plot the geographical map" + "but it is not installed" ) - warnings.warn( - "Basemap will be deprecated in a future release of pysteps, use Cartopy instead", - PendingDeprecationWarning, + # if plot_map == "cartopy": # not really an option for the moment + cartopy_scale = kwargs.get("scale", "50m") + cartopy_subplot = kwargs.get("subplot", (1, 1, 1)) + crs = utils.proj4_to_cartopy(proj4str) + + ax = plot_map_cartopy( + crs, + extent, + cartopy_scale, + drawlonlatlines=drawlonlatlines, + drawlonlatlabels=drawlonlatlabels, + lw=lw, + subplot=cartopy_subplot, ) - ax = Basemap(**bm_params) - - if coastlinecolor is not None: - ax.drawcoastlines(color=coastlinecolor, linewidth=lw, zorder=0.1) - if countrycolor is not None: - ax.drawcountries(color=countrycolor, linewidth=lw, zorder=0.2) - if rivercolor is not None: - ax.drawrivers(zorder=0.2, color=rivercolor) - if continentcolor is not None: - ax.fillcontinents(color=continentcolor, lake_color=lakecolor, zorder=0) - if mapboundarycolor is not None: - ax.drawmapboundary(fill_color=mapboundarycolor, zorder=-1) - if drawlonlatlines: - ax.drawmeridians( - np.linspace(ax.llcrnrlon, ax.urcrnrlon, 10), - color=(0.5, 0.5, 0.5), - linewidth=0.25, - labels=[0, 0, 0, 1], - fmt="%.1f", - fontsize=6, - ) - ax.drawparallels( - np.linspace(ax.llcrnrlat, ax.urcrnrlat, 10), - color=(0.5, 0.5, 0.5), - linewidth=0.25, - labels=[1, 0, 0, 0], - fmt="%.1f", - fontsize=6, - ) - return ax def plot_map_cartopy( - crs, extent, scale, drawlonlatlines=False, lw=0.5, subplot=(1, 1, 1) + crs, + extent, + scale, + drawlonlatlines=False, + drawlonlatlabels=True, + lw=0.5, + subplot=(1, 1, 1), ): """ Plot coastlines, countries, rivers and meridians/parallels using Cartopy. Parameters ---------- - crs : object + crs: object Instance of a crs class defined in cartopy.crs. It can be created using utils.proj4_to_cartopy. - extent : scalars (left, right, bottom, top) + extent: scalars (left, right, bottom, top) The coordinates of the bounding box. - drawlonlatlines : bool + drawlonlatlines: bool Whether to plot longitudes and latitudes. - scale : {'10m', '50m', '110m'} + drawlonlatlabels: bool, optional + If set to True, draw longitude and latitude labels. Valid only if + 'drawlonlatlines' is True. + scale: {'10m', '50m', '110m'} The scale (resolution) of the map. The available options are '10m', '50m', and '110m'. - lw : float + lw: float Line width. - subplot : scalars (nrows, ncols, index) + subplot: scalars (nrows, ncols, index) Subplot dimensions (n_rows, n_cols) and subplot number (index). Returns ------- - ax : axes + ax: axes Cartopy axes. Compatible with matplotlib. """ - if not cartopy_imported: + if not CARTOPY_IMPORTED: raise MissingOptionalDependency( - "plot_map='cartopy' option passed to plot_map_cartopy function " - "but the cartopy package is not installed" + "the cartopy package is required to plot the geographical map" + " but it is not installed" ) if isinstance(subplot, gridspec.SubplotSpec): @@ -333,9 +230,37 @@ def plot_map_cartopy( ), zorder=2, ) + if scale in ["10m", "50m"]: + ax.add_feature( + cfeature.NaturalEarthFeature( + "physical", + "reefs", + scale="10m", + edgecolor="black", + facecolor="none", + linewidth=lw, + ), + zorder=2, + ) + ax.add_feature( + cfeature.NaturalEarthFeature( + "physical", + "minor_islands", + scale="10m", + edgecolor="black", + facecolor="none", + linewidth=lw, + ), + zorder=2, + ) if drawlonlatlines: - ax.gridlines(crs=ccrs.PlateCarree()) + gl = ax.gridlines( + crs=ccrs.PlateCarree(), draw_labels=drawlonlatlabels, dms=True + ) + gl.top_labels = gl.right_labels = False + gl.y_inline = gl.x_inline = False + gl.rotate_labels = False ax.set_extent(extent, crs) diff --git a/pysteps/visualization/motionfields.py b/pysteps/visualization/motionfields.py index 9616e7f4d..b9ef033ac 100644 --- a/pysteps/visualization/motionfields.py +++ b/pysteps/visualization/motionfields.py @@ -13,14 +13,14 @@ import matplotlib.pylab as plt import numpy as np -from pysteps.exceptions import UnsupportedSomercProjection +from pysteps.exceptions import MissingOptionalDependency, UnsupportedSomercProjection from . import basemaps from . import utils def quiver( - UV, ax=None, geodata=None, axis="on", step=20, quiver_kwargs=None, **kwargs, + UV, ax=None, geodata=None, axis="on", step=20, quiver_kwargs={}, map_kwargs={}, ): """Function to plot a motion field as arrows. @@ -34,7 +34,7 @@ def quiver( Array of shape (2,m,n) containing the input motion field. ax : axis object Optional axis object to use for plotting. - geodata : dictionary + geodata : dictionary or None Optional dictionary containing geographical information about the field. @@ -74,7 +74,7 @@ def quiver( Other parameters ---------------- - kwargs: dict + map_kwargs: dict Optional parameters that need to be passed to :py:func:`pysteps.visualization.basemaps.plot_geography`. @@ -84,12 +84,6 @@ def quiver( Figure axes. Needed if one wants to add e.g. text inside the plot. """ - plot_map = kwargs.get("plot_map", None) - if plot_map is not None and geodata is None: - raise ValueError("plot_map!=None but geodata=None") - - if quiver_kwargs is None: - quiver_kwargs = dict() # prepare x y coordinates reproject = False @@ -109,7 +103,7 @@ def quiver( extent = (geodata["x1"], geodata["x2"], geodata["y1"], geodata["y2"]) # check geodata and project if different from axes - if ax is not None and plot_map is None: + if ax is not None: if type(ax).__name__ == "GeoAxesSubplot": try: ccrs = utils.proj4_to_cartopy(geodata["projection"]) @@ -118,8 +112,6 @@ def quiver( # This will work reasonably well for Europe only. t_proj4str = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs" reproject = True - elif type(ax).__name__ == "Basemap": - utils.proj4_to_basemap(geodata["projection"]) if reproject: geodata = utils.reproject_geodata( @@ -136,9 +128,14 @@ def quiver( X, Y = np.meshgrid(x, y) # draw basemaps - if plot_map is not None: + if geodata is not None: try: - ax = basemaps.plot_geography(geodata["projection"], extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], extent, **map_kwargs,) + + except MissingOptionalDependency as e: + # Cartopy is not installed + print(f"{e.__class__}: {e}") + ax = plt.axes() except UnsupportedSomercProjection: # Define default fall-back projection for Swiss data(EPSG:3035) @@ -148,7 +145,8 @@ def quiver( extent = (geodata["x1"], geodata["x2"], geodata["y1"], geodata["y2"]) X, Y = geodata["X_grid"], geodata["Y_grid"] - ax = basemaps.plot_geography(geodata["projection"], extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], extent, **map_kwargs,) + else: ax = plt.gca() @@ -178,7 +176,7 @@ def quiver( def streamplot( - UV, ax=None, geodata=None, axis="on", streamplot_kwargs=None, **kwargs, + UV, ax=None, geodata=None, axis="on", streamplot_kwargs={}, map_kwargs={}, ): """Function to plot a motion field as streamlines. @@ -195,7 +193,7 @@ def streamplot( Array of shape (2, m,n) containing the input motion field. ax : axis object Optional axis object to use for plotting. - geodata : dictionary + geodata : dictionary or None Optional dictionary containing geographical information about the field. If geodata is not None, it must contain the following key-value pairs: @@ -231,7 +229,7 @@ def streamplot( Other parameters ---------------- - kwargs: dict + map_kwargs: dict Optional parameters that need to be passed to :py:func:`pysteps.visualization.basemaps.plot_geography`. @@ -241,13 +239,6 @@ def streamplot( Figure axes. Needed if one wants to add e.g. text inside the plot. """ - plot_map = kwargs.get("plot_map", None) - - if plot_map is not None and geodata is None: - raise ValueError("plot_map!=None but geodata=None") - - if streamplot_kwargs is None: - streamplot_kwargs = dict() # prepare x y coordinates reproject = False @@ -267,7 +258,7 @@ def streamplot( extent = (geodata["x1"], geodata["x2"], geodata["y1"], geodata["y2"]) # check geodata and project if different from axes - if ax is not None and plot_map is None: + if ax is not None: if type(ax).__name__ == "GeoAxesSubplot": try: ccrs = utils.proj4_to_cartopy(geodata["projection"]) @@ -276,8 +267,6 @@ def streamplot( # This will work reasonably well for Europe only. t_proj4str = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs" reproject = True - elif type(ax).__name__ == "Basemap": - utils.proj4_to_basemap(geodata["projection"]) if reproject: geodata = utils.reproject_geodata( @@ -292,9 +281,15 @@ def streamplot( y = np.arange(UV.shape[1]) # draw basemaps - if plot_map is not None: + if geodata is not None: try: - ax = basemaps.plot_geography(geodata["projection"], extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], extent, **map_kwargs,) + + except MissingOptionalDependency as e: + # Cartopy is not installed + print(f"{e.__class__}: {e}") + ax = plt.axes() + except UnsupportedSomercProjection: # Define default fall-back projection for Swiss data(EPSG:3035) # This will work reasonably well for Europe only. @@ -305,7 +300,8 @@ def streamplot( x = X[0, :] y = Y[:, 0] - ax = basemaps.plot_geography(geodata["projection"], extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], extent, **map_kwargs,) + else: ax = plt.gca() diff --git a/pysteps/visualization/precipfields.py b/pysteps/visualization/precipfields.py index 72b515f92..5b45ff239 100755 --- a/pysteps/visualization/precipfields.py +++ b/pysteps/visualization/precipfields.py @@ -39,7 +39,7 @@ def plot_precip_field( colorbar=True, axis="on", cax=None, - **kwargs, + map_kwargs={}, ): """ Function to plot a precipitation intensity or probability field with a @@ -58,7 +58,7 @@ def plot_precip_field( Type of the map to plot: 'intensity' = precipitation intensity field, 'depth' = precipitation depth (accumulation) field, 'prob' = exceedance probability field. - geodata : dictionary, optional + geodata : dictionary or None, optional Optional dictionary containing geographical information about the field. Required is map is not None. @@ -93,9 +93,9 @@ def plot_precip_field( bbox : tuple, optional Four-element tuple specifying the coordinates of the bounding box. Use this for plotting a subdomain inside the input grid. The coordinates are - of the form (lower left x,lower left y,upper right x,upper right y). If - map is not None, the x- and y-coordinates are longitudes and latitudes. - Otherwise they represent image pixels. + of the form (lower left x, lower left y ,upper right x, upper right y). + If 'geodata' is not None, the bbox is in map coordinates, otherwise + it represents image pixels. colorscale : {'pysteps', 'STEPS-BE', 'BOM-RF3'}, optional Which colorscale to use. Applicable if units is 'mm/h', 'mm' or 'dBZ'. probthr : float, optional @@ -114,7 +114,7 @@ def plot_precip_field( Other parameters ---------------- - kwargs: dict + map_kwargs: dict Optional parameters that need to be passed to :py:func:`pysteps.visualization.basemaps.plot_geography`. @@ -135,9 +135,6 @@ def plot_precip_field( ) if type == "prob" and colorbar and probthr is None: raise ValueError("type='prob' but probthr not specified") - plot_map = kwargs.get("plot_map", None) - if plot_map is not None and geodata is None: - raise ValueError("map!=None but geodata=None") if len(R.shape) != 2: raise ValueError("the input is not two-dimensional array") @@ -152,8 +149,8 @@ def plot_precip_field( else: if not PYPROJ_IMPORTED: raise MissingOptionalDependency( - "pyproj package is required to import " - "FMI's radar reflectivity composite " + "pyproj package is required to plot " + "georeferenced precipitation fields " "but it is not installed" ) pr = pyproj.Proj(geodata["projection"]) @@ -166,9 +163,14 @@ def plot_precip_field( origin = "upper" # plot geography - if plot_map is not None: + if geodata is not None: try: - ax = basemaps.plot_geography(geodata["projection"], bm_extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], bm_extent, **map_kwargs,) + regular_grid = True + except MissingOptionalDependency as e: + # Cartopy is not installed + print(f"{e.__class__}: {e}") + ax = plt.axes() regular_grid = True except UnsupportedSomercProjection: # Define default fall-back projection for Swiss data(EPSG:3035) @@ -183,24 +185,21 @@ def plot_precip_field( X, Y = geodata["X_grid"], geodata["Y_grid"] regular_grid = geodata["regular_grid"] - ax = basemaps.plot_geography(geodata["projection"], bm_extent, **kwargs,) + ax = basemaps.plot_geography(geodata["projection"], bm_extent, **map_kwargs,) + else: + ax = plt.axes() regular_grid = True - if bbox is not None and plot_map is not None: + if bbox is not None and geodata is not None: x1, y1 = pr(geodata["x1"], geodata["y1"], inverse=True) x2, y2 = pr(geodata["x2"], geodata["y2"], inverse=True) - if plot_map == "basemap": - x1, y1 = ax(x1, y1) - x2, y2 = ax(x2, y2) - else: - x1, y1 = pr(x1, y1) - x2, y2 = pr(x2, y2) + x1, y1 = pr(x1, y1) + x2, y2 = pr(x2, y2) field_extent = (x1, x2, y1, y2) # plot rainfield if regular_grid: - ax = plt.gca() im = _plot_field( R, ax, type, units, colorscale, extent=field_extent, origin=origin ) @@ -250,19 +249,17 @@ def plot_precip_field( else: cbar.set_label("P(R > %.1f %s)" % (probthr, units)) - if plot_map is None and bbox is not None: - ax = plt.gca() + if geodata is None or axis == "off": + ax.xaxis.set_ticks([]) + ax.xaxis.set_ticklabels([]) + ax.yaxis.set_ticks([]) + ax.yaxis.set_ticklabels([]) + + if bbox is not None: ax.set_xlim(bbox[0], bbox[2]) ax.set_ylim(bbox[1], bbox[3]) - if geodata is None or axis == "off": - axes = plt.gca() - axes.xaxis.set_ticks([]) - axes.xaxis.set_ticklabels([]) - axes.yaxis.set_ticks([]) - axes.yaxis.set_ticklabels([]) - - return plt.gca() + return ax def _plot_field(R, ax, type, units, colorscale, extent, origin=None): @@ -281,16 +278,12 @@ def _plot_field(R, ax, type, units, colorscale, extent, origin=None): else: R[R < 1e-3] = np.nan - vmin, vmax = [None, None] if type in ["intensity", "depth"] else [0.0, 1.0] - im = ax.imshow( R, cmap=cmap, norm=norm, extent=extent, interpolation="nearest", - vmin=vmin, - vmax=vmax, origin=origin, zorder=1, ) diff --git a/pysteps/visualization/utils.py b/pysteps/visualization/utils.py index 2a0a3f204..9277b9f42 100644 --- a/pysteps/visualization/utils.py +++ b/pysteps/visualization/utils.py @@ -8,7 +8,6 @@ :toctree: ../generated/ parse_proj4_string - proj4_to_basemap proj4_to_cartopy reproject_geodata """ @@ -19,15 +18,15 @@ try: import cartopy.crs as ccrs - cartopy_imported = True + CARTOPY_IMPORTED = True except ImportError: - cartopy_imported = False + CARTOPY_IMPORTED = False try: import pyproj - pyproj_imported = True + PYPROJ_IMPORTED = True except ImportError: - pyproj_imported = False + PYPROJ_IMPORTED = False def parse_proj4_string(proj4str): @@ -56,51 +55,6 @@ def parse_proj4_string(proj4str): return result -def proj4_to_basemap(proj4str): - """Convert a PROJ.4 projection string into a dictionary that can be expanded - as keyword arguments to mpl_toolkits.basemap.Basemap.__init__. - - Parameters - ---------- - proj4str : str - A PROJ.4-compatible projection string. - - Returns - ------- - out : dict - The output dictionary. - - """ - pdict = parse_proj4_string(proj4str) - odict = {} - - for k, v in list(pdict.items()): - if k == "proj": - # TODO: Make sure that the proj.4 projection type is in all cases - # mapped to the corresponding (or closest matching) Basemap - # projection. - if v == "somerc": - raise UnsupportedSomercProjection("unsupported projection:" " somerc") - if v not in ["latlon", "latlong", "lonlat", "longlat"]: - odict["projection"] = v - else: - odict["projection"] = "cyl" - elif k == "lon_0" or k == "lat_0" or k == "lat_ts": - # TODO: Check that east/west and north/south hemispheres are - # handled correctly. - if v[-1] in ["E", "N", "S", "W"]: - v = v[:-1] - odict[k] = float(v) - elif k == "ellps": - odict[k] = v - elif k == "R": - odict["rsphere"] = float(v) - elif k in ["k", "k0"]: - odict["k_0"] = float(v) - - return odict - - def proj4_to_cartopy(proj4str): """Convert a PROJ.4 projection string into a Cartopy coordinate reference system (crs) object. @@ -116,13 +70,13 @@ def proj4_to_cartopy(proj4str): Instance of a crs class defined in cartopy.crs. """ - if not cartopy_imported: + if not CARTOPY_IMPORTED: raise MissingOptionalDependency( "cartopy package is required for proj4_to_cartopy function " "utility but it is not installed" ) - if not pyproj_imported: + if not PYPROJ_IMPORTED: raise MissingOptionalDependency( "pyproj package is required for proj4_to_cartopy function utility " "but it is not installed" @@ -182,11 +136,15 @@ def proj4_to_cartopy(proj4str): cl = ccrs.Stereographic elif v == "aea": cl = ccrs.AlbersEqualArea + elif v == "aeqd": + cl = ccrs.AzimuthalEquidistant elif v == "somerc": # Note: ccrs.epsg(2056) doesn't work because the projection # limits are too strict. # We'll use the Stereographic projection as an alternative. cl = ccrs.Stereographic + elif v == "geos": + cl = ccrs.Geostationary else: raise ValueError("unsupported projection: %s" % v) elif k in km_proj: @@ -237,7 +195,7 @@ def reproject_geodata(geodata, t_proj4str, return_grid=None): regular_grid=False to indicate that the reprojected grid has no regular spacing. """ - if not pyproj_imported: + if not PYPROJ_IMPORTED: raise MissingOptionalDependency( "pyproj package is required for reproject_geodata function utility" " but it is not installed"