diff --git a/pysteps/utils/spectral.py b/pysteps/utils/spectral.py
index a304b4abc..a144429b0 100644
--- a/pysteps/utils/spectral.py
+++ b/pysteps/utils/spectral.py
@@ -98,37 +98,43 @@ def mean(X, shape):
 
 
 def rapsd(
-    field, fft_method=None, return_freq=False, d=1.0, normalize=False, **fft_kwargs
+    field, fft_method=np.fft, return_freq=False, d=1.0, normalize=False, **fft_kwargs
 ):
     """
-    Compute radially averaged power spectral density (RAPSD) from the given
-    2D input field.
+    Compute radially averaged power spectral density (RAPSD) of the given 2D
+    input field.
 
     Parameters
     ----------
-    field: array_like
+    field : array_like
         A 2d array of shape (m, n) containing the input field.
-    fft_method: object
+    fft_method : object, optional
+        The fast Fourier transform (FFT) method to apply to the input field.
         A module or object implementing the same methods as numpy.fft and
         scipy.fftpack. If set to None, field is assumed to represent the
         shifted discrete Fourier transform of the input field, where the
         origin is at the center of the array
-        (see numpy.fft.fftshift or scipy.fftpack.fftshift).
-    return_freq: bool
-        Whether to also return the Fourier frequencies.
-    d: scalar
-        Sample spacing (inverse of the sampling rate). Defaults to 1.
-        Applicable if return_freq is 'True'.
-    normalize: bool
-        If True, normalize the power spectrum so that it sums to one.
+        (see numpy.fft.fftshift or scipy.fftpack.fftshift). Defaults to
+        numpy.fft.
+    return_freq : bool, optional
+        Whether to also return the Fourier frequencies. Defaults to ``False``.
+    d : scalar, optional
+        Sample spacing (inverse of the sampling rate). Applicable if
+        return_freq is ``True``. Defaults to 1.
+    normalize : bool, optional
+        If True, normalize the power spectrum so that it sums to one. Defaults
+        to ``False``.
+    fft_kwargs : dict, optional
+        Optional keyword arguments that are passed to the FFT method. Defaults
+        to ``None``.
 
     Returns
     -------
     out: ndarray
-      One-dimensional array containing the RAPSD. The length of the array is
-      int(l/2) (if l is even) or int(l/2)+1 (if l is odd), where l=max(m,n).
+        One-dimensional array containing the RAPSD. The length of the array is
+        int(l/2) (if l is even) or int(l/2)+1 (if l is odd), where l=max(m,n).
     freq: ndarray
-      One-dimensional array containing the Fourier frequencies.
+        One-dimensional array containing the Fourier frequencies.
 
     References
     ----------
diff --git a/pysteps/verification/detcatscores.py b/pysteps/verification/detcatscores.py
index 83a8d1512..3b4201720 100644
--- a/pysteps/verification/detcatscores.py
+++ b/pysteps/verification/detcatscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.detcatscores
 =================================
@@ -27,14 +26,14 @@ def det_cat_fct(pred, obs, thr, scores="", axis=None):
 
     Parameters
     ----------
-    pred: array_like
+    pred : array_like
         Array of predictions. NaNs are ignored.
-    obs: array_like
+    obs : array_like
         Array of verifying observations. NaNs are ignored.
-    thr: float
+    thr : float
         The threshold that is applied to predictions and observations in order
-        to define events vs no events (yes/no).
-    scores: {string, list of strings}, optional
+        to define events vs. no events (yes/no).
+    scores : {string, list of strings}, optional
         The name(s) of the scores. The default, scores="", will compute all
         available scores.
         The available score names are:
@@ -73,7 +72,7 @@ def det_cat_fct(pred, obs, thr, scores="", axis=None):
         |  SEDI      | symmetric extremal dependency index                    |
         +------------+--------------------------------------------------------+
 
-    axis: None or int or tuple of ints, optional
+    axis : None or int or tuple of ints, optional
         Axis or axes along which a score is integrated. The default, axis=None,
         will integrate all of the elements of the input arrays.\n
         If axis is -1 (or any negative integer),
@@ -84,16 +83,17 @@ def det_cat_fct(pred, obs, thr, scores="", axis=None):
 
     Returns
     -------
-    result: dict
+    result : dict
         Dictionary containing the verification results.
 
     See also
     --------
-    pysteps.verification.detcontscores.det_cont_fct
+    :py:func:`pysteps.verification.detcatscores.det_cont_fct`
     """
 
     contab = det_cat_fct_init(thr, axis)
     det_cat_fct_accum(contab, pred, obs)
+
     return det_cat_fct_compute(contab, scores)
 
 
@@ -103,10 +103,10 @@ def det_cat_fct_init(thr, axis=None):
 
     Parameters
     ----------
-    thr: float
+    thr : float
         threshold that is applied to predictions and observations in order
-        to define events vs no events (yes/no).
-    axis: None or int or tuple of ints, optional
+        to define events vs. no events (yes/no).
+    axis : None or int or tuple of ints, optional
         Axis or axes along which a score is integrated. The default, axis=None,
         will integrate all of the elements of the input arrays.\n
         If axis is -1 (or any negative integer),
@@ -118,7 +118,7 @@ def det_cat_fct_init(thr, axis=None):
     Returns
     -------
     out: dict
-      The contingency table object.
+        The contingency table object.
     """
 
     contab = {}
@@ -148,12 +148,12 @@ def det_cat_fct_accum(contab, pred, obs):
 
     Parameters
     ----------
-    contab: dict
-      A contingency table object initialized with
-      pysteps.verification.detcatscores.det_cat_fct_init.
-    pred: array_like
+    contab : dict
+        A contingency table object initialized with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`.
+    pred : array_like
         Array of predictions. NaNs are ignored.
-    obs: array_like
+    obs : array_like
         Array of verifying observations. NaNs are ignored.
     """
 
@@ -222,21 +222,21 @@ def det_cat_fct_merge(contab_1, contab_2):
 
     Parameters
     ----------
-    contab_1: dict
-      A contingency table object initialized with
-      :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`
-      and populated with
-      :py:func:`pysteps.verification.detcatscores.det_cat_fct_accum`.
-    contab_2: dict
-      Another contingency table object initialized with
-      :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`
-      and populated with
-      :py:func:`pysteps.verification.detcatscores.det_cat_fct_accum`.
+    contab_1 : dict
+        A contingency table object initialized with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`
+        and populated with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_accum`.
+    contab_2 : dict
+        Another contingency table object initialized with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`
+        and populated with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_accum`.
 
     Returns
     -------
-    out: dict
-      The merged contingency table object.
+    out : dict
+        The merged contingency table object.
     """
 
     # checks
@@ -270,11 +270,12 @@ def det_cat_fct_compute(contab, scores=""):
 
     Parameters
     ----------
-    contab: dict
-      A contingency table object initialized with
-      pysteps.verification.detcatscores.det_cat_fct_init and populated with
-      pysteps.verification.detcatscores.det_cat_fct_accum.
-    scores: {string, list of strings}, optional
+    contab : dict
+        A contingency table object initialized with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_init`
+        and populated with
+        :py:func:`pysteps.verification.detcatscores.det_cat_fct_accum`.
+    scores : {string, list of strings}, optional
         The name(s) of the scores. The default, scores="", will compute all
         available scores.
         The available score names a
@@ -315,7 +316,7 @@ def det_cat_fct_compute(contab, scores=""):
 
     Returns
     -------
-    result: dict
+    result : dict
         Dictionary containing the verification results.
     """
 
diff --git a/pysteps/verification/detcontscores.py b/pysteps/verification/detcontscores.py
index 25d989016..d34eb7545 100644
--- a/pysteps/verification/detcontscores.py
+++ b/pysteps/verification/detcontscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.detcontscores
 ==================================
@@ -26,11 +25,11 @@ def det_cont_fct(pred, obs, scores="", axis=None, conditioning=None, thr=0.0):
 
     Parameters
     ----------
-    pred: array_like
+    pred : array_like
         Array of predictions. NaNs are ignored.
-    obs: array_like
+    obs : array_like
         Array of verifying observations. NaNs are ignored.
-    scores: {string, list of strings}, optional
+    scores : {string, list of strings}, optional
         The name(s) of the scores. The default, scores="", will compute all
         available scores.
         The available score names are:
@@ -69,7 +68,7 @@ def det_cont_fct(pred, obs, scores="", axis=None, conditioning=None, thr=0.0):
         |            | as in Germann et al. (2006)                            |
         +------------+--------------------------------------------------------+
 
-    axis: {int, tuple of int, None}, optional
+    axis : {int, tuple of int, None}, optional
         Axis or axes along which a score is integrated. The default, axis=None,
         will integrate all of the elements of the input arrays.\n
         If axis is -1 (or any negative integer),
@@ -77,18 +76,18 @@ def det_cont_fct(pred, obs, scores="", axis=None, conditioning=None, thr=0.0):
         and scores are computed on all of the elements in the input arrays.\n
         If axis is a tuple of ints, the integration is performed on all of the
         axes specified in the tuple.
-    conditioning: {None, "single", "double"}, optional
+    conditioning : {None, "single", "double"}, optional
         The type of conditioning used for the verification.
         The default, conditioning=None, includes all pairs. With
         conditioning="single", only pairs with either pred or obs > thr are
         included. With conditioning="double", only pairs with both pred and
         obs > thr are included.
-    thr: float
+    thr : float
         Optional threshold value for conditioning. Defaults to 0.
 
     Returns
     -------
-    result: dict
+    result : dict
         Dictionary containing the verification results.
 
     Notes
@@ -126,7 +125,7 @@ def det_cont_fct(pred, obs, scores="", axis=None, conditioning=None, thr=0.0):
 
     See also
     --------
-    pysteps.verification.detcatscores.det_cat_fct
+    :py:func:`pysteps.verification.detcatscores.det_cat_fct`
     """
 
     # catch case of single score passed as string
@@ -211,7 +210,7 @@ def det_cont_fct_init(axis=None, conditioning=None, thr=0.0):
 
     Parameters
     ----------
-    axis: {int, tuple of int, None}, optional
+    axis : {int, tuple of int, None}, optional
         Axis or axes along which a score is integrated. The default, axis=None,
         will integrate all of the elements of the input arrays.\n
         If axis is -1 (or any negative integer),
@@ -219,18 +218,18 @@ def det_cont_fct_init(axis=None, conditioning=None, thr=0.0):
         and scores are computed on all of the elements in the input arrays.\n
         If axis is a tuple of ints, the integration is performed on all of the
         axes specified in the tuple.
-    conditioning: {None, "single", "double"}, optional
+    conditioning : {None, "single", "double"}, optional
         The type of conditioning used for the verification.
         The default, conditioning=None, includes all pairs. With
         conditioning="single", only pairs with either pred or obs > thr are
         included. With conditioning="double", only pairs with both pred and
         obs > thr are included.
-    thr: float
+    thr : float, optional
         Optional threshold value for conditioning. Defaults to 0.
 
     Returns
     -------
-    out: dict
+    out : dict
         The verification error object.
     """
 
@@ -267,12 +266,12 @@ def det_cont_fct_accum(err, pred, obs):
 
     Parameters
     ----------
-    err: dict
+    err : dict
         A verification error object initialized with
         :py:func:`pysteps.verification.detcontscores.det_cont_fct_init`.
-    pred: array_like
+    pred : array_like
         Array of predictions. NaNs are ignored.
-    obs: array_like
+    obs : array_like
         Array of verifying observations. NaNs are ignored.
 
     References
@@ -396,21 +395,21 @@ def det_cont_fct_merge(err_1, err_2):
 
     Parameters
     ----------
-    err_1: dict
-      A verification error object initialized with
-      :py:func:`pysteps.verification.detcontscores.det_cont_fct_init`
-      and populated with
-      :py:func:`pysteps.verification.detcontscores.det_cont_fct_accum`.
-    err_2: dict
-      Another verification error object initialized with
-      :py:func:`pysteps.verification.detcontscores.det_cont_fct_init`
-      and populated with
-      :py:func:`pysteps.verification.detcontscores.det_cont_fct_accum`.
+    err_1 : dict
+        A verification error object initialized with
+        :py:func:`pysteps.verification.detcontscores.det_cont_fct_init`
+        and populated with
+        :py:func:`pysteps.verification.detcontscores.det_cont_fct_accum`.
+    err_2 : dict
+        Another verification error object initialized with
+        :py:func:`pysteps.verification.detcontscores.det_cont_fct_init`
+        and populated with
+        :py:func:`pysteps.verification.detcontscores.det_cont_fct_accum`.
 
     Returns
     -------
-    out: dict
-      The merged verification error object.
+    out : dict
+        The merged verification error object.
     """
 
     # checks
@@ -522,7 +521,7 @@ def det_cont_fct_compute(err, scores=""):
 
     Returns
     -------
-    result: dict
+    result : dict
         Dictionary containing the verification results.
     """
 
diff --git a/pysteps/verification/ensscores.py b/pysteps/verification/ensscores.py
index df6b00088..b32a3d6f0 100644
--- a/pysteps/verification/ensscores.py
+++ b/pysteps/verification/ensscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.ensscores
 ==============================
diff --git a/pysteps/verification/interface.py b/pysteps/verification/interface.py
index 6e325f6f4..dd3d054a8 100644
--- a/pysteps/verification/interface.py
+++ b/pysteps/verification/interface.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.interface
 ==============================
diff --git a/pysteps/verification/lifetime.py b/pysteps/verification/lifetime.py
index 4c5dcb03f..cd8f0c501 100644
--- a/pysteps/verification/lifetime.py
+++ b/pysteps/verification/lifetime.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.lifetime
 =============================
diff --git a/pysteps/verification/plots.py b/pysteps/verification/plots.py
index 21974b013..d726b653a 100755
--- a/pysteps/verification/plots.py
+++ b/pysteps/verification/plots.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.plots
 ==========================
diff --git a/pysteps/verification/probscores.py b/pysteps/verification/probscores.py
index f1eb398ba..b8ded142d 100644
--- a/pysteps/verification/probscores.py
+++ b/pysteps/verification/probscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.probscores
 ===============================
@@ -31,17 +30,17 @@ def CRPS(X_f, X_o):
 
     Parameters
     ----------
-    X_f: array_like
-      Array of shape (k,m,n,...) containing the values from an ensemble
-      forecast of k members with shape (m,n,...).
-    X_o: array_like
-      Array of shape (m,n,...) containing the observed values corresponding
-      to the forecast.
+    X_f : array_like
+        Array of shape (k,m,n,...) containing the values from an ensemble
+        forecast of k members with shape (m,n,...).
+    X_o : array_like
+        Array of shape (m,n,...) containing the observed values corresponding
+        to the forecast.
 
     Returns
     -------
-    out: float
-      The computed CRPS.
+    out : float
+        The computed CRPS.
 
     References
     ----------
@@ -61,8 +60,8 @@ def CRPS_init():
 
     Returns
     -------
-    out: dict
-      The CRPS object.
+    out : dict
+        The CRPS object.
     """
     return {"CRPS_sum": 0.0, "n": 0.0}
 
@@ -75,14 +74,15 @@ def CRPS_accum(CRPS, X_f, X_o):
 
     Parameters
     ----------
-    CRPS: dict
-      The CRPS object.
-    X_f: array_like
-      Array of shape (k,m,n,...) containing the values from an ensemble
-      forecast of k members with shape (m,n,...).
-    X_o: array_like
-      Array of shape (m,n,...) containing the observed values corresponding
-      to the forecast.
+    CRPS : dict
+        The CRPS object initialized with
+        :py:func:`pysteps.verification.probscores.CRPS_init`.
+    X_f : array_like
+        Array of shape (k,m,n,...) containing the values from an ensemble
+        forecast of k members with shape (m,n,...).
+    X_o : array_like
+        Array of shape (m,n,...) containing the observed values corresponding
+        to the forecast.
 
     References
     ----------
@@ -137,47 +137,50 @@ def CRPS_compute(CRPS):
 
     Parameters
     ----------
-    CRPS: dict
-      A CRPS object created with CRPS_init.
+    CRPS : dict
+        A CRPS object initialized with
+        :py:func:`pysteps.verification.probscores.CRPS_init`.
 
     Returns
     -------
-    out: float
-      The computed CRPS.
+    out : float
+        The computed CRPS.
     """
     return 1.0 * CRPS["CRPS_sum"] / CRPS["n"]
 
 
 def reldiag(P_f, X_o, X_min, n_bins=10, min_count=10):
     """
-    Compute the x- and y- coordinates of the points in the reliability diagram.
+    Compute a reliability diagram from the given forecast exceedance
+    probabilities and observations.
 
     Parameters
     ----------
-    P_f: array-like
-      Forecast probabilities for exceeding the intensity threshold specified
-      in the reliability diagram object.
-    X_o: array-like
-      Observed values.
-    X_min: float
-      Precipitation intensity threshold for yes/no prediction.
-    n_bins: int
-        Number of bins to use in the reliability diagram.
-    min_count: int
-      Minimum number of samples required for each bin. A zero value is assigned
-      if the number of samples in a bin is smaller than bin_count.
+    P_f : array_like
+        Forecast probabilities for exceeding the intensity threshold ``X_min``.
+    X_o : array_like
+        Observed values that are converted into yes/no predictions of exceeding
+        the threshold ``X_min``.
+    X_min : float
+        Precipitation intensity threshold for yes/no prediction.
+    n_bins: int, optional
+        Number of bins to use in the reliability diagram. Defaults to 10.
+    min_count: int, optional
+        Minimum number of samples required for each bin. A zero value is
+        assigned if the number of samples in a bin is smaller than
+        ``bin_count``. Defaults to 10.
 
     Returns
     -------
-    out: tuple
-      Two-element tuple containing the x- and y-coordinates of the points in
-      the reliability diagram.
+    out : tuple
+        Two-element tuple containing the x- and y-coordinates of the points in
+        the reliability curve.
     """
-
     P_f = P_f.copy()
     X_o = X_o.copy()
     rdiag = reldiag_init(X_min, n_bins, min_count)
     reldiag_accum(rdiag, P_f, X_o)
+
     return reldiag_compute(rdiag)
 
 
@@ -187,18 +190,19 @@ def reldiag_init(X_min, n_bins=10, min_count=10):
 
     Parameters
     ----------
-    X_min: float
-      Precipitation intensity threshold for yes/no prediction.
-    n_bins: int
-        Number of bins to use in the reliability diagram.
-    min_count: int
-      Minimum number of samples required for each bin. A zero value is assigned
-      if the number of samples in a bin is smaller than bin_count.
+    X_min : float
+        Precipitation intensity threshold for yes/no prediction.
+    n_bins : int, optional
+        Number of bins to use in the reliability diagram. Defaults to 10.
+    min_count : int, optional
+        Minimum number of samples required for each bin. A zero value is
+        assigned if the number of samples in a bin is smaller than
+        ``bin_count``. Defaults to 10.
 
     Returns
     -------
-    out: dict
-      The reliability diagram object.
+    out : dict
+        The reliability diagram object.
 
     References
     ----------
@@ -219,18 +223,20 @@ def reldiag_init(X_min, n_bins=10, min_count=10):
 
 
 def reldiag_accum(reldiag, P_f, X_o):
-    """Accumulate the given probability-observation pairs into the reliability
+    """Accumulate the given probability-observation pairs into a reliability
     diagram.
 
     Parameters
     ----------
-    reldiag: dict
-      A reliability diagram object created with reldiag_init.
-    P_f: array-like
-      Forecast probabilities for exceeding the intensity threshold specified
-      in the reliability diagram object.
-    X_o: array-like
-      Observed values.
+    reldiag : dict
+        A reliability diagram object initialized with
+        :py:func:`pysteps.verification.probscores.reldiag_init`.
+    P_f : array-like
+        Forecast probabilities for exceeding the intensity threshold specified
+        in the reliability diagram object.
+    X_o : array-like
+        Observed values that are converted into yes/no predictions of exceeding
+        the intensity threshold specified in the reliability diagram object.
     """
     mask = np.logical_and(np.isfinite(P_f), np.isfinite(X_o))
 
@@ -266,18 +272,20 @@ def reldiag_accum(reldiag, P_f, X_o):
 
 def reldiag_compute(reldiag):
     """
-    Compute the x- and y- coordinates of the points in the reliability diagram.
+    Compute the x- and y- coordinates of the curve points from the given
+    reliability diagram object.
 
     Parameters
     ----------
-    reldiag: dict
-      A reliability diagram object created with reldiag_init.
+    reldiag : dict
+        A reliability diagram object initialized with
+        :py:func:`pysteps.verification.probscores.reldiag_init`.
 
     Returns
     -------
-    out: tuple
-      Two-element tuple containing the x- and y-coordinates of the points in
-      the reliability diagram.
+    out : tuple
+        Two-element tuple containing the x- and y-coordinates of the points in
+        the reliability curve.
     """
     f = 1.0 * reldiag["Y_sum"] / reldiag["num_idx"]
     r = 1.0 * reldiag["X_sum"] / reldiag["num_idx"]
@@ -287,36 +295,42 @@ def reldiag_compute(reldiag):
 
 def ROC_curve(P_f, X_o, X_min, n_prob_thrs=10, compute_area=False):
     """
-    Compute the ROC curve and its area from the given ROC object.
+    Compute the ROC curve from the given forecast exceedance probabilities and
+    observations.
 
     Parameters
     ----------
-    P_f: array_like
-      Forecasted probabilities for exceeding the threshold specified in the ROC
-      object. Non-finite values are ignored.
-    X_o: array_like
-      Observed values. Non-finite values are ignored.
-    X_min: float
-      Precipitation intensity threshold for yes/no prediction.
-    n_prob_thrs: int
-      The number of probability thresholds to use.
-      The interval [0,1] is divided into n_prob_thrs evenly spaced values.
-    compute_area: bool
-      If True, compute the area under the ROC curve (between 0.5 and 1).
+    P_f : array_like
+        Forecasted probabilities for exceeding the threshold ``X_min``.
+        Non-finite values are ignored.
+    X_o : array_like
+        Observed values that are compared against the threshold ``X_min`` and
+        converted into yes/no predictions of exceedance. Non-finite values are
+        ignored.
+    X_min : float
+        Precipitation intensity threshold for yes/no prediction of exceedance.
+    n_prob_thrs : int, optional
+        The number of probability thresholds to use. The interval [0,1] is
+        divided into ``n_prob_thrs`` evenly spaced values. Defaults to 10.
+    compute_area : bool, optional
+        If ``True``, compute the area under the ROC curve (between 0.5 and 1).
+        Defaults to ``False``.
 
     Returns
     -------
-    out: tuple
-      A two-element tuple containing the probability of detection (POD) and
-      probability of false detection (POFD) for the probability thresholds
-      specified in the ROC curve object. If compute_area is True, return the
-      area under the ROC curve as the third element of the tuple.
+    out : tuple
+        A two-element tuple containing the probabilities of false detection
+        (POFD, x-coordinate) and the probabilities of detection (POD,
+        y-coordinate) for the specified probability thresholds. If compute_area
+        is ``True``, return the area under the ROC curve as the third element
+        of the tuple.
     """
 
     P_f = P_f.copy()
     X_o = X_o.copy()
     roc = ROC_curve_init(X_min, n_prob_thrs)
     ROC_curve_accum(roc, P_f, X_o)
+
     return ROC_curve_compute(roc, compute_area)
 
 
@@ -326,16 +340,16 @@ def ROC_curve_init(X_min, n_prob_thrs=10):
 
     Parameters
     ----------
-    X_min: float
-      Precipitation intensity threshold for yes/no prediction.
-    n_prob_thrs: int
-      The number of probability thresholds to use.
-      The interval [0,1] is divided into n_prob_thrs evenly spaced values.
+    X_min : float
+        Precipitation intensity threshold for yes/no prediction of exceedance.
+    n_prob_thrs : int, optional
+        The number of probability thresholds to use. The interval [0,1] is
+        divided into ``n_prob_thrs`` evenly spaced values. Defaults to 10.
 
     Returns
     -------
-    out: dict
-      The ROC curve object.
+    out : dict
+        The ROC curve object.
     """
     ROC = {}
 
@@ -351,17 +365,20 @@ def ROC_curve_init(X_min, n_prob_thrs=10):
 
 def ROC_curve_accum(ROC, P_f, X_o):
     """Accumulate the given probability-observation pairs into the given ROC
-    object.
+    curve object.
 
     Parameters
     ----------
-    ROC: dict
-      A ROC curve object created with ROC_curve_init.
-    P_f: array_like
-      Forecasted probabilities for exceeding the threshold specified in the ROC
-      object. Non-finite values are ignored.
-    X_o: array_like
-      Observed values. Non-finite values are ignored.
+    ROC : dict
+        A ROC curve object initialized with
+        :py:func:`pysteps.verification.probscores.ROC_curve_init`.
+    P_f : array_like
+        Forecasted probabilities for exceeding the threshold specified in the
+        ROC object. Non-finite values are ignored.
+    X_o : array_like
+        Observed values that are compared against the intensity threshold
+        specified in the ROC curve object and converted into yes/no predictions
+        of exceedance. Non-finite values are ignored.
     """
     mask = np.logical_and(np.isfinite(P_f), np.isfinite(X_o))
 
@@ -385,18 +402,21 @@ def ROC_curve_compute(ROC, compute_area=False):
 
     Parameters
     ----------
-    ROC: dict
-      A ROC curve object created with ROC_curve_init.
-    compute_area: bool
-      If True, compute the area under the ROC curve (between 0.5 and 1).
+    ROC : dict
+        A ROC curve object initialized with
+        :py:func:`pysteps.verification.probscores.ROC_curve_init`.
+    compute_area : bool, optional
+        If ``True``, compute the area under the ROC curve (between 0.5 and 1).
+        Defaults to ``False``.
 
     Returns
     -------
-    out: tuple
-      A two-element tuple containing the probability of detection (POD) and
-      probability of false detection (POFD) for the probability thresholds
-      specified in the ROC curve object. If compute_area is True, return the
-      area under the ROC curve as the third element of the tuple.
+    out : tuple
+        A two-element tuple containing the probabilities of false detection
+        (POFD, x-coordinate) and the probabilities of detection (POD,
+        y-coordinate) for the probability thresholds specified in the ROC curve
+        object. If compute_area is ``True``, return the area under the ROC
+        curve as the third element of the tuple.
     """
     POD_vals = []
     POFD_vals = []
diff --git a/pysteps/verification/salscores.py b/pysteps/verification/salscores.py
index f2d97b77c..1f1f04f07 100644
--- a/pysteps/verification/salscores.py
+++ b/pysteps/verification/salscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.salscores
 ==============================
diff --git a/pysteps/verification/spatialscores.py b/pysteps/verification/spatialscores.py
index 41e02c93e..e0f48d316 100644
--- a/pysteps/verification/spatialscores.py
+++ b/pysteps/verification/spatialscores.py
@@ -1,4 +1,3 @@
-# -- coding: utf-8 --
 """
 pysteps.verification.spatialscores
 ==================================