Set new LK defaults (#151)

* Set default value for epsilon from 5 to 10
* Set default value for buffer_mask from 0 to 5
* Change function name to lowercase to comply with pep8
* Adapt test for interfaces
* Clean up docstrings

Author: dnerini

examples/LK_buffer_mask.py

@@ -97,22 +97,22 @@
 # Sparse Lucas-Kanade
 # -------------------
 #
-# By setting the optional argument 'dense=False' in 'xy, uv = LK_optflow(.....)',
+# By setting the optional argument ``dense=False`` in ``xy, uv = dense_lucaskanade(...)``,
 # the LK algorithm returns the motion vectors detected by the Lucas-Kanade scheme
 # without interpolating them on the grid.
 # This allows us to better identify the presence of wrongly detected
 # stationary motion in areas where precipitation is leaving the domain (look
 # for the red dots within the blue circle in the figure below).
 
-# get Lucas-Kanade optical flow method
-LK_optflow = motion.get_method("LK")
+# Get Lucas-Kanade optical flow method
+dense_lucaskanade = motion.get_method("LK")
 
 # Mask invalid values
 R = np.ma.masked_invalid(R)
 
-# Use default settings (i.e., no buffering of the radar mask)
-fd_kwargs1 = {"buffer_mask":0}
-xy, uv = LK_optflow(R, dense=False, fd_kwargs=fd_kwargs1)
+# Use no buffering of the radar mask
+fd_kwargs1 = {"buffer_mask" : 0}
+xy, uv = dense_lucaskanade(R, dense=False, fd_kwargs=fd_kwargs1)
 plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys"))
 plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5)
 plt.quiver(
@@ -127,22 +127,23 @@
 )
 circle = plt.Circle((620, 245), 100, color="b", clip_on=False, fill=False)
 plt.gca().add_artist(circle)
-plt.title("buffer_mask = 0 (default)")
+plt.title("buffer_mask = 0")
 plt.show()
 
 ################################################################################
-# By default, the LK algorithm considers missing values as no precipitation, i.e.,
+# The LK algorithm cannot distinguish missing values from no precipitation, that is,
 # no-data are the same as no-echoes. As a result, the fixed boundaries produced
 # by precipitation in contact with no-data areas are interpreted as stationary motion.
 # One way to mitigate this effect of the boundaries is to introduce a slight buffer
 # of the no-data mask so that the algorithm will ignore all the portions of the
 # radar domain that are nearby no-data areas.
-# This is achieved by passing the keyword argument 'buffer_mask = 10' within the
-# feature detection optional arguments 'fd_kwargs'.
+# This buffer can be set by the keyword argument ``buffer_mask`` within the
+# feature detection optional arguments ``fd_kwargs``.
+# Note that by default ``dense_lucaskanade`` uses a 5-pixel buffer.
 
 # with buffer
-buffer = 10
-fd_kwargs2 = {"buffer_mask":buffer}
+buffer = 5
+fd_kwargs2 = {"buffer_mask" : buffer}
 xy, uv = LK_optflow(R, dense=False, fd_kwargs=fd_kwargs2)
 plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys"))
 plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5)
@@ -166,13 +167,16 @@
 # ------------------
 #
 # The above displacement vectors produced by the Lucas-Kanade method are now
-# interpolated to produce a full field of motion (i.e., 'dense=True').
+# interpolated to produce a full field of motion (i.e., ``dense=True``).
 # By comparing the velocity of the motion fields, we can easily notice
 # the negative bias that is introduced by the the erroneous interpretation of
 # velocities near the maximum range of the radars.
+# Please note that we are setting a small shape parameter ``epsilon`` for the
+# interpolation routine. This will produce a smoother motion field.
 
-UV1 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs1)
-UV2 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs2)
+interp_kwargs = {"epsilon" : 5} # use a small shape parameter for interpolation
+UV1 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs1, interp_kwargs=interp_kwargs)
+UV2 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs2, interp_kwargs=interp_kwargs)
 
 V1 = np.sqrt(UV1[0] ** 2 + UV1[1] ** 2)
 V2 = np.sqrt(UV2[0] ** 2 + UV2[1] ** 2)
@@ -183,16 +187,16 @@
 plt.show()
 
 ################################################################################
-# Notice that the default motion field can be significantly slower (more than 10%
-# slower) because of the presence of wrong stationary motion vectors at the edges.
+# Notice how the presence of erroneous velocity vectors produces a significantly
+# slower motion field near the right edge of the domain.
 #
 # Forecast skill
 # --------------
 #
 # We are now going to evaluate the benefit of buffering the radar mask by computing
 # the forecast skill in terms of the Spearman correlation coefficient.
 # The extrapolation forecasts are computed using the dense UV motion fields
-# estimted above.
+# estimated above.
 
 # Get the advection routine and extrapolate the last radar frame by 12 time steps
 # (i.e., 1 hour lead time)

pysteps/tests/test_interfaces.py

@@ -229,7 +229,7 @@ def test_utils_interface():
                              ('clip', dimension.clip_domain),
                              ('square', dimension.square_domain),
                              ('upscale', dimension.aggregate_fields_space),
-                             ('shitomasi', images.ShiTomasi_detection),
+                             ('shitomasi', images.shitomasi_detection),
                              ('morph_opening', images.morph_opening),
                              ('rbfinterp2d', interpolate.rbfinterp2d),
                              ('rapsd', spectral.rapsd),

pysteps/utils/cleansing.py

@@ -31,12 +31,12 @@ def decluster(coord, input_array, scale, min_samples=1, verbose=False):
     input_array : array_like
         Array of shape (n) or (n, m), where *n* is the number of samples and
         *m* the number of variables.
-        All values in **input_array** are required to have finite values.
+        All values in ``input_array`` are required to have finite values.
 
     scale : float or array_like
-        The **scale** parameter in the same units of **coord**.
+        The ``scale`` parameter in the same units of ``coord``.
         It can be a scalar or an array_like of shape (d).
-        Data points within the declustering **scale** are aggregated.
+        Data points within the declustering ``scale`` are aggregated.
 
     min_samples : int, optional
         The minimum number of samples for computing the median within a given
@@ -49,8 +49,8 @@ def decluster(coord, input_array, scale, min_samples=1, verbose=False):
     -------
 
     out : tuple of ndarrays
-        A two-element tuple (**out_coord**, **output_array**) containing the
-        declustered coordinates (l, d) and **input_array** (l, m), where *l* is
+        A two-element tuple (``out_coord``, ``output_array``) containing the
+        declustered coordinates (l, d) and input array (l, m), where *l* is
         the new number of samples with *l* <= *n*.
     """
 
@@ -147,22 +147,21 @@ def detect_outliers(input_array, thr, coord=None, k=None, verbose=False):
         Array of shape (n) or (n, m), where *n* is the number of samples and
         *m* the number of variables. If *m* > 1, the Mahalanobis distance
         is used.
-        All values in **input_array** are required to have finite values.
+        All values in ``input_array`` are required to have finite values.
 
     thr : float
-        The number of standard deviations from the mean
-        that defines an outlier.
+        The number of standard deviations from the mean used to define an outlier.
 
     coord : array_like, optional
         Array of shape (n, d) containing the coordinates of the input data into
         a space of *d* dimensions.
-        Passing **coord** requires that **k** is not None.
+        Passing ``coord`` requires that ``k`` is not None.
 
     k : int or None, optional
         The number of nearest neighbours used to localize the outlier
         detection.
         If set to None (the default), it employs all the data points (global
-        detection). Setting **k** requires that **coord** is not None.
+        detection). Setting ``k`` requires that ``coord`` is not None.
 
     verbose : bool, optional
         Print out information.
@@ -172,7 +171,7 @@ def detect_outliers(input_array, thr, coord=None, k=None, verbose=False):
 
     out : array_like
         A 1-D boolean array of shape (n) with True values indicating the outliers
-        detected in **input_array**.
+        detected in ``input_array``.
     """
 
     input_array = np.copy(input_array)

pysteps/utils/images.py

@@ -12,7 +12,7 @@
 .. autosummary::
     :toctree: ../generated/
 
-    ShiTomasi_detection
+    shitomasi_detection
     morph_opening
 """
 
@@ -29,12 +29,12 @@
     CV2_IMPORTED = False
 
 
-def ShiTomasi_detection(input_image,
+def shitomasi_detection(input_image,
                         max_corners=1000,
                         quality_level=0.01,
                         min_distance=10,
                         block_size=5,
-                        buffer_mask=0,
+                        buffer_mask=5,
                         use_harris=False,
                         k=0.04,
                         verbose=False,
@@ -75,27 +75,27 @@ def ShiTomasi_detection(input_image,
         valid pixels.
 
     max_corners : int, optional
-        The **maxCorners** parameter in the `Shi-Tomasi`_ corner detection
+        The ``maxCorners`` parameter in the `Shi-Tomasi`_ corner detection
         method.
         It represents the maximum number of points to be tracked (corners).
         If set to zero, all detected corners are used.
 
     quality_level : float, optional
-        The **qualityLevel** parameter in the `Shi-Tomasi`_ corner detection
+        The ``qualityLevel`` parameter in the `Shi-Tomasi`_ corner detection
         method.
         It represents the minimal accepted quality for the image corners.
 
     min_distance : int, optional
-        The **minDistance** parameter in the `Shi-Tomasi`_ corner detection
+        The ``minDistance`` parameter in the `Shi-Tomasi`_ corner detection
         method.
         It represents minimum possible Euclidean distance in pixels between
         corners.
 
     block_size : int, optional
-        The **blockSize** parameter in the `Shi-Tomasi`_ corner detection
+        The ``blockSize`` parameter in the `Shi-Tomasi`_ corner detection
         method.
         It represents the window size in pixels used for computing a derivative
-        covariation matrix over each pixel neighborhood.
+        covariation matrix over each pixel neighbourhood.
 
     use_harris : bool, optional
         Whether to use a `Harris detector`_  or cornerMinEigenVal_.

pysteps/utils/interface.py

@@ -179,7 +179,7 @@ def donothing(R, metadata=None, *args, **kwargs):
     methods_objects["upscale"] = dimension.aggregate_fields_space
 
     # image processing methods
-    methods_objects["shitomasi"] = images.ShiTomasi_detection
+    methods_objects["shitomasi"] = images.shitomasi_detection
     methods_objects["morph_opening"] = images.morph_opening
 
     # interpolation methods

pysteps/utils/interpolate.py

@@ -22,7 +22,7 @@ def rbfinterp2d(
     xgrid,
     ygrid,
     rbfunction="gaussian",
-    epsilon=5,
+    epsilon=10,
     k=50,
     nchunks=5,
 ):
@@ -43,27 +43,27 @@ def rbfinterp2d(
         Array of shape (n) or (n, m) containing the values of the data points,
         where *n* is the number of data points and *m* the number of co-located
         variables.
-        All values in **input_array** are required to have finite values.
+        All values in ``input_array`` are required to have finite values.
 
     xgrid, ygrid : array_like
         1D arrays representing the coordinates of the 2-D output grid.
 
-    rbfunction : {"gaussian", "multiquadric", "inverse quadratic",
-        "inverse multiquadric", "bump"}, optional
+    rbfunction : {"gaussian", "multiquadric", "inverse quadratic", "inverse multiquadric", "bump"}, optional
         The name of one of the available radial basis function based on a
-        normalized Euclidian norm.
+        normalized Euclidian norm as defined in the **Notes** section below.
 
-        See also the Notes section below.
+        More details provided in the wikipedia reference page linked below.
 
     epsilon : float, optional
         The shape parameter used to scale the input to the radial kernel.
 
-        A smaller value for **epsilon** produces a smoother interpolation. More
-        details provided in the wikipedia reference page.
+        A smaller value for ``epsilon`` produces a smoother interpolation. More
+        details provided in the wikipedia reference page linked below.
 
     k : int or None, optional
-        The number of nearest neighbours used to speed-up the interpolation.
-        If set to None, it interpolates based on all the data points.
+        The number of nearest neighbours used for each target location.
+        This can also be useful to to speed-up the interpolation.
+        If set to None, it interpolates using all the data points at once.
 
     nchunks : int, optional
         The number of chunks in which the grid points are split to limit the
@@ -73,7 +73,7 @@ def rbfinterp2d(
     -------
 
     output_array : ndarray_
-        The interpolated field(s) having shape (m, ygrid.size, xgrid.size).
+        The interpolated field(s) having shape (*m*, ``ygrid.size``, ``xgrid.size``).
 
     Notes
     -----