Refactor Optical flow interfaces (#110)

* Add decorator to check the input_images shape in all optical flow functions
* Set private LK methods as public
* Simplify syntax, improve docstrings
* Remove extra_vectors argument in LK
* Refactor outlier detection
* Refactor declustering
* Refactor interpolation routine
* Use numpy.empty instead of empty lists
* Improve usage of MaskedArray
* Add utils.cleansing, utils.images and utils.interpolate modules
* Update tests
* Update doc

Author: dnerini

doc/source/pysteps_reference/utils.rst

@@ -7,8 +7,11 @@ Implementation of miscellaneous utility functions.
 
 .. automodule:: pysteps.utils.interface
 .. automodule:: pysteps.utils.arrays
+.. automodule:: pysteps.utils.cleansing
 .. automodule:: pysteps.utils.conversion
 .. automodule:: pysteps.utils.dimension
 .. automodule:: pysteps.utils.fft
+.. automodule:: pysteps.utils.images
+.. automodule:: pysteps.utils.interpolate
 .. automodule:: pysteps.utils.spectral
 .. automodule:: pysteps.utils.transformation

examples/LK_buffer_mask.py

@@ -1,15 +1,16 @@
+# -*- coding: utf-8 -*-
 """
 Handling of no-data in Lucas-Kanade
 ===================================
 
 Areas of missing data in radar images are typically caused by visibility limits
-such as beam blockage and the radar coverage itself. These artifacts can mislead 
+such as beam blockage and the radar coverage itself. These artifacts can mislead
 the echo tracking algorithms. For instance, precipitation leaving the domain
 might be erroneously detected as having nearly stationary velocity.
 
-This example shows how the Lucas-Kanade algorithm can be tuned to avoid the 
-erroneous interpretation of velocities near the maximum range of the radars by 
-buffering the no-data mask in the radar image in order to exclude all vectors 
+This example shows how the Lucas-Kanade algorithm can be tuned to avoid the
+erroneous interpretation of velocities near the maximum range of the radars by
+buffering the no-data mask in the radar image in order to exclude all vectors
 detected nearby no-data areas.
 """
 
@@ -71,7 +72,9 @@
 mask[~np.isnan(ref_mm)] = np.nan
 
 # Log-transform the data [dBR]
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+R, metadata = transformation.dB_transform(
+    R, metadata, threshold=0.1, zerovalue=-15.0
+)
 
 # Keep the reference frame in dBR (for plotting purposes)
 ref_dbr = R[0].copy()
@@ -109,10 +112,20 @@
 R.data[R.mask] = np.nan
 
 # Use default settings (i.e., no buffering of the radar mask)
-x, y, u, v = LK_optflow(R, dense=False, buffer_mask=0, quality_level_ST=0.1)
+fd_kwargs1 = {"buffer_mask":0}
+xy, uv = LK_optflow(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(x, y, u, v, color="red", angles="xy", scale_units="xy", scale=0.2)
+plt.quiver(
+    xy[:, 0],
+    xy[:, 1],
+    uv[:, 0],
+    uv[:, 1],
+    color="red",
+    angles="xy",
+    scale_units="xy",
+    scale=0.2,
+)
 circle = plt.Circle((620, 245), 100, color="b", clip_on=False, fill=False)
 plt.gca().add_artist(circle)
 plt.title("buffer_mask = 0 (default)")
@@ -129,13 +142,24 @@
 # 'x,y,u,v = LK_optflow(.....)'.
 
 # with buffer
-x, y, u, v = LK_optflow(R, dense=False, buffer_mask=20, quality_level_ST=0.2)
+buffer = 10
+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)
-plt.quiver(x, y, u, v, color="red", angles="xy", scale_units="xy", scale=0.2)
+plt.quiver(
+    xy[:, 0],
+    xy[:, 1],
+    uv[:, 0],
+    uv[:, 1],
+    color="red",
+    angles="xy",
+    scale_units="xy",
+    scale=0.2,
+)
 circle = plt.Circle((620, 245), 100, color="b", clip_on=False, fill=False)
 plt.gca().add_artist(circle)
-plt.title("buffer_mask = 20")
+plt.title("buffer_mask = %i" % buffer)
 plt.show()
 
 ################################################################################
@@ -148,13 +172,13 @@
 # the negative bias that is introduced by the the erroneous interpretation of
 # velocities near the maximum range of the radars.
 
-UV1 = LK_optflow(R, dense=True, buffer_mask=0, quality_level_ST=0.1)
-UV2 = LK_optflow(R, dense=True, buffer_mask=20, quality_level_ST=0.2)
+UV1 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs1)
+UV2 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs2)
 
 V1 = np.sqrt(UV1[0] ** 2 + UV1[1] ** 2)
 V2 = np.sqrt(UV2[0] ** 2 + UV2[1] ** 2)
 
-plt.imshow((V1 - V2) / V2, cmap=cm.RdBu_r, vmin=-0.1, vmax=0.1)
+plt.imshow((V1 - V2) / V2, cmap=cm.RdBu_r, vmin=-0.5, vmax=0.5)
 plt.colorbar(fraction=0.04, pad=0.04)
 plt.title("Relative difference in motion speed")
 plt.show()
@@ -184,7 +208,13 @@
 
 # Find the veriyfing observations in the archive
 fns = io.archive.find_by_date(
-    date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=12
+    date,
+    root_path,
+    path_fmt,
+    fn_pattern,
+    fn_ext,
+    timestep=5,
+    num_next_files=12,
 )
 
 # Read and convert the radar composites
@@ -200,8 +230,8 @@
     score_2.append(skill(R_f2[i, :, :], R_o[i + 1, :, :])["corr_s"])
 
 x = (np.arange(12) + 1) * 5  # [min]
-plt.plot(x, score_1, label="no mask buffer")
-plt.plot(x, score_2, label="with mask buffer")
+plt.plot(x, score_1, label="buffer_mask = 0")
+plt.plot(x, score_2, label="buffer_mask = %i" % buffer)
 plt.legend()
 plt.xlabel("Lead time [min]")
 plt.ylabel("Corr. coeff. []")

pysteps/decorators.py

@@ -0,0 +1,56 @@
+"""
+pysteps.decorators
+==================
+
+Decorators used to define reusable building blocks that can change or extend
+the behavior of some functions in pysteps.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    check_motion_input_image
+"""
+from functools import wraps
+
+import numpy as np
+
+
+def check_input_frames(minimum_input_frames=2,
+                       maximum_input_frames=np.inf,
+                       just_ndim=False):
+    """
+    Check that the input_images used as inputs in the optical-flow
+    methods has the correct shape (t, x, y ).
+    """
+
+    def _check_input_frames(motion_method_func):
+        @wraps(motion_method_func)
+        def new_function(*args, **kwargs):
+            """
+            Return new function with the checks prepended to the
+            target motion_method_func function.
+            """
+
+            input_images = args[0]
+            if input_images.ndim != 3:
+                raise ValueError(
+                    "input_images dimension mismatch.\n"
+                    f"input_images.shape: {str(input_images.shape)}\n"
+                    "(t, x, y ) dimensions expected"
+                )
+
+            if not just_ndim:
+                num_of_frames = input_images.shape[0]
+
+                if minimum_input_frames < num_of_frames > maximum_input_frames:
+                    raise ValueError(
+                        f"input_images frames {num_of_frames} mismatch.\n"
+                        f"Minimum frames: {minimum_input_frames}\n"
+                        f"Maximum frames: {maximum_input_frames}\n"
+                    )
+
+            return motion_method_func(*args, **kwargs)
+
+        return new_function
+
+    return _check_input_frames