Refactor visualization module (#199)

1) Refactor visualization.utils functions 

- parse_proj4_string: A more robust implementation of this function is implemented in pyproj. The refactored implementation uses a direct call to the pyproj function.
- proj4_to_cartopy: Refactor the proj4_to_cartopy function
  - Use the parse_proj4_string to construct a dictionary with the proj definitions
  - Simplify the function code by using additional dictionaries to relate the Pyproj and Cartopy parameters. 
  - Add support for additional cartopy keywords defined in the Proj4 definitions.
- Add two plotting utility functions to avoid code duplication and to improve the axis handling:
  - get_basemap_axis: This function safely gets a basemap axis. If the current axis is not a cartopy axis, it creates a basemap. Otherwise, the current axis is returned.
  - get_geogrid: Returns the x,y grid used in the plots along with additional metadata. 
-  Solves issue #197.

2) Refactor plotting functions

- Refactor the plotting functions to use the new plotting utils
- Merge `_plot_field_pcolormesh` and `_plot_field` in a single function.
- Avoid plotting the no-data mask separately. Instead, plot it in the `_plot_field` function. 
- Merge `quiver` and `streamplots` functions into a single one called `plot_motion` for avoiding code duplication.
The `quiver` and `streamplot` functions are now wrappers for `plot_motion`.
- **IMPORTANT**
  - Rename the `type` keyword in  `plot_precip_field` to `ptype` to avoid shadowing the built-in function with the same name. Add the corresponding deprecation warning.

3) Define a unique zorder for each type of plot

- Basemap features zorder
  - ocean: 0
  - land: 0
  - lakes: 0
  - rivers_lake_centerlines: 0
  - coastline: 15
  - cultural: 15
  - reefs: 15
  - minor_islands: 15
- precipitation: 10
- quiver: 20
- stream function: 30
- tracks and contours: 40

Co-authored-by: Daniel Nerini

Author: aperezhortal

examples/advection_correction.py

@@ -17,7 +17,7 @@
 """
 
 from datetime import datetime
-import matplotlib.pyplot as pl
+import matplotlib.pyplot as plt
 import numpy as np
 
 from pysteps import io, motion, rcparams
@@ -134,13 +134,13 @@ def advection_correction(R, T=5, t=1):
 # correction to account for this spatial shift. The final result is a smoother
 # rainfall accumulation map.
 
-pl.figure(figsize=(9, 4))
-pl.subplot(121)
+plt.figure(figsize=(9, 4))
+plt.subplot(121)
 plot_precip_field(R.mean(axis=0), title="3-h rainfall accumulation")
-pl.subplot(122)
+plt.subplot(122)
 plot_precip_field(R_ac, title="Same with advection correction")
-pl.tight_layout()
-pl.show()
+plt.tight_layout()
+plt.show()
 
 ################################################################################
 # Reference

examples/plot_cascade_decomposition.py

@@ -8,7 +8,7 @@
 
 """
 
-from matplotlib import cm, pyplot
+from matplotlib import cm, pyplot as plt
 import numpy as np
 import os
 from pprint import pprint
@@ -40,6 +40,7 @@
 
 # Plot the rainfall field
 plot_precip_field(R, geodata=metadata)
+plt.show()
 
 # Log-transform the data
 R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
@@ -58,14 +59,15 @@
 
 # Plot the power spectrum
 M, N = F.shape
-fig, ax = pyplot.subplots()
+fig, ax = plt.subplots()
 im = ax.imshow(
     np.log(F ** 2), vmin=4, vmax=24, cmap=cm.jet, extent=(-N / 2, N / 2, -M / 2, M / 2)
 )
 cb = fig.colorbar(im)
 ax.set_xlabel("Wavenumber $k_x$")
 ax.set_ylabel("Wavenumber $k_y$")
 ax.set_title("Log-power spectrum of R")
+plt.show()
 
 ###############################################################################
 # Cascade decomposition
@@ -81,7 +83,7 @@
 
 # Plot the bandpass filter weights
 L = max(N, M)
-fig, ax = pyplot.subplots()
+fig, ax = plt.subplots()
 for k in range(num_cascade_levels):
     ax.semilogx(
         np.linspace(0, L / 2, len(filter["weights_1d"][k, :])),
@@ -97,6 +99,7 @@
 ax.set_xlabel("Radial wavenumber $|\mathbf{k}|$")
 ax.set_ylabel("Normalized weight")
 ax.set_title("Bandpass filter weights")
+plt.show()
 
 ###############################################################################
 # Finally, apply the 2D Gaussian filters to decompose the radar rainfall field
@@ -110,7 +113,7 @@
     sigma = decomp["stds"][i]
     decomp["cascade_levels"][i] = (decomp["cascade_levels"][i] - mu) / sigma
 
-fig, ax = pyplot.subplots(nrows=2, ncols=4)
+fig, ax = plt.subplots(nrows=2, ncols=4)
 
 ax[0, 0].imshow(R, cmap=cm.RdBu_r, vmin=-5, vmax=5)
 ax[0, 1].imshow(decomp["cascade_levels"][0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
@@ -134,6 +137,7 @@
     for j in range(4):
         ax[i, j].set_xticks([])
         ax[i, j].set_yticks([])
-pyplot.tight_layout()
+plt.tight_layout()
+plt.show()
 
 # sphinx_gallery_thumbnail_number = 4

examples/plot_noise_generators.py

@@ -11,7 +11,7 @@
 field.
 """
 
-from matplotlib import cm, pyplot
+from matplotlib import cm, pyplot as plt
 import numpy as np
 import os
 from pprint import pprint
@@ -44,6 +44,7 @@
 
 # Plot the rainfall field
 plot_precip_field(R, geodata=metadata)
+plt.show()
 
 # Log-transform the data
 R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
@@ -79,7 +80,7 @@
 b2 = Fp["pars"][3]
 
 # Plot the observed power spectrum and the model
-fig, ax = pyplot.subplots()
+fig, ax = plt.subplots()
 plot_scales = [512, 256, 128, 64, 32, 16, 8, 4]
 plot_spectrum1d(
     freq,
@@ -101,11 +102,12 @@
     label="Fit",
     wavelength_ticks=plot_scales,
 )
-pyplot.legend()
+plt.legend()
 ax.set_title(
     "Radially averaged log-power spectrum of R\n"
     r"$\omega_0=%.0f km, \beta_1=%.1f, \beta_2=%.1f$" % (w0, b1, b2)
 )
+plt.show()
 
 ###############################################################################
 # Nonparametric filter
@@ -136,7 +138,7 @@
 
 # Plot the generated noise fields
 
-fig, ax = pyplot.subplots(nrows=2, ncols=3)
+fig, ax = plt.subplots(nrows=2, ncols=3)
 
 # parametric noise
 ax[0, 0].imshow(Np[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
@@ -152,7 +154,8 @@
     for j in range(3):
         ax[i, j].set_xticks([])
         ax[i, j].set_yticks([])
-pyplot.tight_layout()
+plt.tight_layout()
+plt.show()
 
 ###############################################################################
 # The above figure highlights the main limitation of the parametric approach

examples/plot_steps_nowcast.py

@@ -61,6 +61,7 @@
 
 # Plot the rainfall field
 plot_precip_field(R[-1, :, :], geodata=metadata)
+plt.show()
 
 # Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
 # set the fill value to -15 dBR
@@ -104,6 +105,7 @@
     geodata=metadata,
     title="S-PROG (+ %i min)" % (n_leadtimes * timestep),
 )
+plt.show()
 
 ###############################################################################
 # As we can see from the figure above, the forecast produced by S-PROG is a
@@ -150,6 +152,7 @@
     geodata=metadata,
     title="Ensemble mean (+ %i min)" % (n_leadtimes * timestep),
 )
+plt.show()
 
 ###############################################################################
 # The mean of the ensemble displays similar properties as the S-PROG
@@ -166,6 +169,7 @@
     )
     ax.set_title("Member %02d" % i)
 plt.tight_layout()
+plt.show()
 
 ###############################################################################
 # As we can see from these two members of the ensemble, the stochastic forecast
@@ -183,10 +187,11 @@
 plot_precip_field(
     P,
     geodata=metadata,
-    type="prob",
+    ptype="prob",
     units="mm/h",
     probthr=0.5,
     title="Exceedence probability (+ %i min)" % (n_leadtimes * timestep),
 )
+plt.show()
 
 # sphinx_gallery_thumbnail_number = 5

examples/rainfarm_downscale.py

@@ -46,8 +46,8 @@
 pprint(metadata)
 
 # Plot the original rainfall field
-plt.figure()
 plot_precip_field(precip, geodata=metadata)
+plt.show()
 
 # Assign the fill value to all the Nans
 precip[~np.isfinite(precip)] = metadata["zerovalue"]
@@ -116,7 +116,7 @@
     plt.subplots_adjust(wspace=0, hspace=0)
 
 plt.tight_layout()
-
+plt.show()
 
 ###############################################################################
 # Remarks

pysteps/exceptions.py

@@ -19,14 +19,3 @@ class DataModelError(Exception):
     """Raised when a file is not compilant with the Data Information Model."""
 
     pass
-
-
-class UnsupportedSomercProjection(Exception):
-    """
-    Raised when the Swiss Oblique Mercator (somerc) projection is passed
-    to cartopy.
-    Necessary since cartopy doesn't support the Swiss projection.
-    TODO: remove once the somerc projection is supported in cartopy.
-    """
-
-    pass

pysteps/tests/test_plt_animate.py

@@ -0,0 +1,54 @@
+# -*- coding: utf-8 -*-
+
+import os
+
+import numpy as np
+import pytest
+
+from pysteps.tests.helpers import get_precipitation_fields
+from pysteps.visualization.animations import animate
+
+
+def test_animate(tmp_path):
+
+    # test data
+    precip, metadata = get_precipitation_fields(
+        num_prev_files=2,
+        num_next_files=0,
+        return_raw=True,
+        metadata=True,
+        upscale=2000,
+    )
+
+    # obs only
+    animate(precip)
+    animate(precip, title="title")
+    animate(precip, timestamps_obs=metadata["timestamps"])
+    animate(precip, geodata=metadata, map_kwargs={"plot_map": None})
+    with pytest.raises(ValueError):
+        animate(precip, timestamps_obs=metadata["timestamps"][:2])
+    animate(precip, motion_field=np.ones((2, *precip.shape[1:])))
+    with pytest.raises(ValueError):
+        animate(precip, motion_plot="test")
+
+    # with forecast
+    animate(precip, precip)
+    animate(precip, precip, title="title")
+    animate(precip, precip, timestamps_obs=metadata["timestamps"])
+    animate(precip, precip, timestamps_obs=metadata["timestamps"], timestep_min=5)
+    with pytest.raises(ValueError):
+        animate(precip, precip, ptype="prob")
+    animate(precip, precip, ptype="prob", prob_thr=1)
+    animate(precip, precip, ptype="mean")
+    animate(precip, np.stack((precip, precip)), ptype="ensemble")
+
+    # save frames
+    animate(
+        precip,
+        np.stack((precip, precip)),
+        display_animation=False,
+        savefig=True,
+        path_outputs=tmp_path,
+        fig_dpi=10,
+    )
+    assert len(os.listdir(tmp_path)) == 9

pysteps/tests/test_plt_cartopy.py

@@ -32,7 +32,7 @@ def test_visualization_plot_precip_field(source, map_kwargs, pass_geodata):
     if not pass_geodata:
         metadata = None
 
-    plot_precip_field(field, type="intensity", geodata=metadata, map_kwargs=map_kwargs)
+    plot_precip_field(field, ptype="intensity", geodata=metadata, map_kwargs=map_kwargs)
 
 
 if __name__ == "__main__":

pysteps/tests/test_plt_precipfields.py

@@ -60,7 +60,7 @@ def test_visualization_plot_precip_field(
     field_orig = field.copy()
     ax = plot_precip_field(
         field.copy(),
-        type=plot_type,
+        ptype=plot_type,
         bbox=bbox,
         geodata=None,
         colorscale=colorscale,