Repository: fatiando/verde
Pull request #58 (2018-05-22 09:12:29)

Make block reduce into a filter class
The new `BlockReduce` class implements the `filter` method and can be
used in a `Chain`.


Patch:
diff --git a/examples/block_reduce.py b/examples/block_reduce.py
index 316a661..60e3da1 100644
--- a/examples/block_reduce.py
+++ b/examples/block_reduce.py
@@ -4,7 +4,7 @@ Data decimation in blocks
 
 When gridding data that has been highly oversampled in a direction (shipborne
 and airborne data, for example), it is important to decimate the data before
-interpolation to avoid aliasing. Function :func:`verde.block_reduce` decimates
+interpolation to avoid aliasing. Class :func:`verde.BlockReduce` decimates
 data by applying a reduction operation (mean, median, mode, max, etc) to the
 data in blocks. For non-smooth data, like bathymetry, a blocked median filter
 is a good choice.
@@ -19,9 +19,9 @@ import verde as vd
 data = vd.datasets.fetch_baja_bathymetry()
 
 # Decimate the data using a blocked median with 10 arc-minute blocks
-coordinates, bathymetry = vd.block_reduce(
-    coordinates=(data.longitude, data.latitude), data=data.bathymetry_m,
-    reduction=np.median, spacing=10/60)
+reducer = vd.BlockReduce(reduction=np.median, spacing=10/60)
+coordinates, bathymetry = reducer.filter((data.longitude, data.latitude),
+                                         data.bathymetry_m)
 lon, lat = coordinates
 
 print("Original data size:", data.bathymetry_m.size)
diff --git a/examples/chain_trend.py b/examples/chain_trend.py
index 415436d..bed07a6 100644
--- a/examples/chain_trend.py
+++ b/examples/chain_trend.py
@@ -3,11 +3,18 @@ Chaining operations
 ===================
 
 The :class:`verde.Chain` class allows us to created gridders that perform
-multiple estimator operations on data. It will fit each estimator on the data
-residuals of the previous one. For example, say we want to calculate a
-polynomial trend for our data, fit a gridder on the residual values (not the
-trend), but then make a grid of the original data values. This is useful
-because many gridders can't handle trends in data.
+multiple operations on data. Each step in the chain filters the input and
+passes the result along to the next step. For gridders and trend estimators,
+filtering means fitting the model and passing along the residuals (input data
+minus predicted data). When predicting data, the predictions of each step are
+added together. Other operations, like :class:`verde.BlockReduce` change
+the input data values and the coordinates as well but don't impact the
+predictions.
+
+For example, say we want to decimate some data, calculate a polynomial trend,
+and fit a gridder on the residual values (not the trend), but then make a grid
+of the full data values (trend included). This is useful because decimation is
+common and many gridders can't handle trends in data.
 """
 import numpy as np
 import matplotlib.pyplot as plt
@@ -21,25 +28,23 @@ import verde as vd
 data = vd.datasets.fetch_rio_magnetic_anomaly()
 region = vd.get_region((data.longitude, data.latitude))
 
-# Before gridding, we need to decimate the data using a blocked mean to avoid
-# aliasing
-spacing = 0.5/60
-coordinates, total_field = vd.block_reduce(
-    (data.longitude, data.latitude), data.total_field_anomaly_nt,
-    reduction=np.median, spacing=spacing)
-
-# Project the data using pyproj so that we can use it as input for the gridder.
-# We'll set the latitude of true scale to the mean latitude of the data.
+# Create a projection for the data using pyproj so that we can use it as input
+# for the gridder. We'll set the latitude of true scale to the mean latitude of
+# the data.
 projection = pyproj.Proj(proj='merc', lat_ts=data.latitude.mean())
-proj_coordinates = projection(*coordinates)
 
-# Create a chain that fits a 2nd degree trend to the anomaly data and then a
-# standard gridder to the residuals
-chain = vd.Chain([('trend', vd.Trend(degree=2)),
-                  ('spline', vd.Spline(damping=1e-8))])
+# Create a chain that decimates the data using a blocked mean to avoid
+# aliasing, fits a 2nd degree trend, and then fits a standard gridder to the
+# residuals. The spacing for the blocked mean will be 0.5 arc-minutes
+# (approximately converted to meters)
+spacing = 0.5/60
+chain = vd.Chain([('reduce', vd.BlockReduce(np.mean, spacing*111e3)),
+                  ('trend', vd.Trend(degree=2)),
+                  ('spline', vd.Spline(damping=1e-7))])
 print("Chained estimator:", chain)
 # Calling 'fit' will automatically run the data through the chain
-chain.fit(proj_coordinates, total_field)
+chain.fit(projection(data.longitude.values, data.latitude.values),
+          data.total_field_anomaly_nt)
 
 # Each component of the chain can be accessed separately using the
 # 'named_steps' attribute
@@ -65,15 +70,13 @@ crs = ccrs.PlateCarree()
 
 plt.figure(figsize=(7, 5))
 ax = plt.axes(projection=ccrs.Mercator())
-ax.set_title("Chained trend and spline")
+ax.set_title("Decimate, de-trend, and spline")
 maxabs = np.max(np.abs([grid.total_field_anomaly.min(),
                         grid.total_field_anomaly.max()]))
 pc = ax.pcolormesh(grid.longitude, grid.latitude, grid.total_field_anomaly,
                    transform=crs, cmap='seismic', vmin=-maxabs, vmax=maxabs)
 cb = plt.colorbar(pc, pad=0.01)
 cb.set_label('total field anomaly [nT]')
-# Plot the locations of the decimated data
-ax.plot(*coordinates, '.k', markersize=0.5, transform=crs)
 # Set the proper ticks for a Cartopy map
 ax.set_xticks(np.arange(-42.5, -42, 0.1), crs=crs)
 ax.set_yticks(np.arange(-22.4, -22, 0.1), crs=crs)
diff --git a/examples/scipygridder.py b/examples/scipygridder.py
index a285608..df68a10 100644
--- a/examples/scipygridder.py
+++ b/examples/scipygridder.py
@@ -31,9 +31,9 @@ data = vd.datasets.fetch_baja_bathymetry()
 # the oversampling along the ship tracks. We'll use a blocked median with 5
 # arc-minute blocks.
 spacing = 5/60
-coordinates, bathymetry = vd.block_reduce(
-    (data.longitude, data.latitude), data.bathymetry_m, reduction=np.median,
-    spacing=spacing)
+reducer = vd.BlockReduce(reduction=np.median, spacing=spacing)
+coordinates, bathymetry = reducer.filter((data.longitude, data.latitude),
+                                         data.bathymetry_m)
 
 # Project the data using pyproj so that we can use it as input for the gridder.
 # We'll set the latitude of true scale to the mean latitude of the data.
diff --git a/examples/spline.py b/examples/spline.py
index f75d0fc..9c27083 100644
--- a/examples/spline.py
+++ b/examples/spline.py
@@ -25,9 +25,9 @@ region = vd.get_region((data.longitude, data.latitude))
 
 # Before gridding, we need to decimate the data to avoid aliasing.
 spacing = 10/60
-coordinates, bathymetry = vd.block_reduce(
-    (data.longitude, data.latitude), data.bathymetry_m, reduction=np.median,
-    spacing=spacing)
+reducer = vd.BlockReduce(reduction=np.median, spacing=spacing)
+coordinates, bathymetry = reducer.filter((data.longitude, data.latitude),
+                                         data.bathymetry_m)
 
 # Now we can add some outliers to test our spline
 outliers = np.random.RandomState(2).randint(0, bathymetry.size, size=20)
diff --git a/verde/__init__.py b/verde/__init__.py
index 7318024..400814e 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -8,7 +8,7 @@ from . import datasets
 from .base import BaseGridder
 from .coordinates import scatter_points, grid_coordinates, inside, \
     profile_coordinates, get_region
-from .grid_math import block_reduce, block_split, distance_mask
+from .grid_math import block_split, distance_mask, BlockReduce
 from .scipy_bridge import ScipyGridder
 from .trend import Trend, trend_jacobian, VectorTrend
 from .chain import Chain
diff --git a/verde/grid_math.py b/verde/grid_math.py
index 39a469f..2abba1e 100644
--- a/verde/grid_math.py
+++ b/verde/grid_math.py
@@ -4,6 +4,7 @@ Operations on spatial data: block operations, derivatives, etc.
 import numpy as np
 import pandas as pd
 from scipy.spatial import cKDTree  # pylint: disable=no-name-in-module
+from sklearn.base import BaseEstimator
 
 from .coordinates import get_region, grid_coordinates
 
@@ -45,7 +46,7 @@ def block_split(coordinates, spacing, adjust='spacing', region=None):
 
     See also
     --------
-    block_reduce : Apply a reduction operation to the data in blocks (windows).
+    BlockReduce : Apply a reduction operation to the data in blocks (windows).
 
     Examples
     --------
@@ -79,19 +80,17 @@ def block_split(coordinates, spacing, adjust='spacing', region=None):
     return block_coords, labels
 
 
-def block_reduce(coordinates, data, reduction, spacing, region=None,
-                 adjust='spacing', center_coordinates=False):
+class BlockReduce(BaseEstimator):
     """
-    Apply a reduction operation to the data in blocks (windows).
+    Apply a reduction/aggregation operation to the data in blocks (windows).
 
     Returns the reduced data value for each block along with the associated
     coordinates, which can be determined through the same reduction applied to
     the coordinates or as the center of each block.
 
     If a data region to be divided into blocks is not given, it will be the
-    bounding region of the data. When using this function to decimate data
-    before gridding, it's best to use the same region and spacing as the
-    desired grid.
+    bounding region of the data. When using this class to decimate data before
+    gridding, it's best to use the same region and spacing as the desired grid.
 
     If the given region is not divisible by the spacing (block size), either
     the region or the spacing will have to be adjusted. By default, the spacing
@@ -100,14 +99,12 @@ def block_reduce(coordinates, data, reduction, spacing, region=None,
 
     Blocks without any data are omitted from the output.
 
+    Implements the :meth:`~verde.BlockReduce.filter` method so it can be used
+    with :class:`verde.Chain`. Only acts during data fitting and is ignored
+    during prediction.
+
     Parameters
     ----------
-    coordinates : tuple of arrays
-        Arrays with the coordinates of each data point. Should be in the
-        following order: (easting, northing, vertical, ...). Only easting and
-        northing will be used, all subsequent coordinates will be ignored.
-    data : array
-        The data values at each point.
     reduction : function
         A reduction function that takes an array and returns a single value
         (e.g., ``np.mean``, ``np.median``, etc).
@@ -127,35 +124,70 @@ def block_reduce(coordinates, data, reduction, spacing, region=None,
         block. Otherwise, the coordinates are calculated by applying the same
         reduction operation to the input coordinates.
 
-    Returns
-    -------
-    blocked_coordinates : tuple of arrays
-        (easting, northing) arrays with the coordinates of each block that
-        contains data.
-    blocked_data : array
-        The block reduced data values.
-
     See also
     --------
     block_split : Split a region into blocks and label points accordingly.
+    verde.Chain : Apply filter operations successively on data.
 
     """
-    easting, northing = coordinates[:2]
-    block_coords, labels = block_split((easting, northing), spacing, adjust,
-                                       region)
-    if center_coordinates:
-        table = pd.DataFrame(dict(data=data.ravel(), block=labels))
-        blocked = table.groupby('block').aggregate(reduction)
-        unique = np.unique(labels)
-        blocked_coords = tuple(i[unique] for i in block_coords)
-    else:
-        table = pd.DataFrame(dict(easting=easting.ravel(),
-                                  northing=northing.ravel(),
-                                  data=data.ravel(),
-                                  block=labels))
-        blocked = table.groupby('block').aggregate(reduction)
-        blocked_coords = (blocked.easting.values, blocked.northing.values)
-    return blocked_coords, blocked.data.values
+
+    def __init__(self, reduction, spacing, region=None, adjust='spacing',
+                 center_coordinates=False):
+        self.reduction = reduction
+        self.spacing = spacing
+        self.region = region
+        self.adjust = adjust
+        self.center_coordinates = center_coordinates
+
+    def filter(self, coordinates, data, weights=None):
+        """
+        Apply the blocked aggregation to the given data.
+
+        Returns the reduced data value for each block along with the associated
+        coordinates, which can be determined through the same reduction applied
+        to the coordinates or as the center of each block.
+
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting
+            and northing will be used, all subsequent coordinates will be
+            ignored.
+        data : array
+            The data values at each point.
+        weights : None or array or tuple of arrays
+            If not None, then the weights assigned to each data point. If more
+            than one data component is provided, you must provide a weights
+            array for each data component (if not None).
+
+        Returns
+        -------
+        blocked_coordinates : tuple of arrays
+            (easting, northing) arrays with the coordinates of each block that
+            contains data.
+        blocked_data : array
+            The block reduced data values.
+
+        """
+        if weights is not None:
+            raise NotImplementedError()
+        easting, northing = coordinates[:2]
+        block_coords, labels = block_split((easting, northing), self.spacing,
+                                           self.adjust, self.region)
+        if self.center_coordinates:
+            table = pd.DataFrame(dict(data=data.ravel(), block=labels))
+            blocked = table.groupby('block').aggregate(self.reduction)
+            unique = np.unique(labels)
+            blocked_coords = tuple(i[unique] for i in block_coords)
+        else:
+            table = pd.DataFrame(dict(easting=easting.ravel(),
+                                      northing=northing.ravel(),
+                                      data=data.ravel(),
+                                      block=labels))
+            blocked = table.groupby('block').aggregate(self.reduction)
+            blocked_coords = (blocked.easting.values, blocked.northing.values)
+        return blocked_coords, blocked.data.values
 
 
 def distance_mask(coordinates, data_coordinates, maxdist):


Test patch:
diff --git a/verde/tests/test_grid_math.py b/verde/tests/test_grid_math.py
index 1d43066..c6d46f0 100644
--- a/verde/tests/test_grid_math.py
+++ b/verde/tests/test_grid_math.py
@@ -3,9 +3,10 @@ Test the grid math functions
 """
 import numpy as np
 import numpy.testing as npt
+import pytest
 
 from ..coordinates import grid_coordinates, scatter_points
-from ..grid_math import block_reduce
+from ..grid_math import BlockReduce
 
 
 def test_block_reduce():
@@ -13,8 +14,10 @@ def test_block_reduce():
     region = (-5, 0, 5, 10)
     east, north = grid_coordinates(region, spacing=0.1, pixel_register=True)
     data = 20*np.ones_like(east)
-    block_coords, block_data = block_reduce((east, north), data, np.mean,
-                                            spacing=1)
+    reducer = BlockReduce(np.mean, spacing=1)
+    with pytest.raises(NotImplementedError):
+        reducer.filter((east, north), data, weights=np.ones_like(data))
+    block_coords, block_data = reducer.filter((east, north), data)
     assert len(block_coords[0]) == 25
     assert len(block_coords[1]) == 25
     assert len(block_data) == 25
@@ -28,9 +31,9 @@ def test_block_reduce_scatter():
     region = (-5, 0, 5, 10)
     east, north = scatter_points(region, size=10000, random_state=0)
     data = 20*np.ones_like(east)
-    block_coords, block_data = block_reduce((east, north), data, np.mean,
-                                            spacing=1, region=region,
-                                            center_coordinates=True)
+    block_coords, block_data = BlockReduce(
+        np.mean, 1, region=region, center_coordinates=True).filter(
+            (east, north), data)
     assert len(block_coords[0]) == 25
     assert len(block_coords[1]) == 25
     assert len(block_data) == 25