Repository: fatiando/verde
Pull request #184 (2019-05-14 01:10:52)

Allow block_split and BlockReduce to take a shape
`block_split` and `BlockReduce`/`BlockMean` (which relied on it) only
allowed defining the block size by `spacing`. Now they an take a `shape`
as well. Useful when the size of the block is less meaningful than the
number of blocks (like in spatial cross-validation splitting).



**Reminders**

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst` and `verde/__init__.py`.
- [ ] Write detailed docstrings for all functions/classes/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.


Patch:
diff --git a/verde/blockreduce.py b/verde/blockreduce.py
index 687b09e..2ececfe 100644
--- a/verde/blockreduce.py
+++ b/verde/blockreduce.py
@@ -39,6 +39,10 @@ class BlockReduce(BaseEstimator):  # pylint: disable=too-few-public-methods
     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.
 
+    The size of the blocks can be specified by the *spacing* parameter. Alternatively,
+    the number of blocks in the South-North and West-East directions can be specified
+    using the *shape* parameter.
+
     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
     will be rounded to the nearest multiple. Optionally, the East and North
@@ -55,6 +59,8 @@ class BlockReduce(BaseEstimator):  # pylint: disable=too-few-public-methods
     reduction : function
         A reduction function that takes an array and returns a single value
         (e.g., ``np.mean``, ``np.median``, etc).
+    shape : tuple = (n_north, n_east) or None
+        The number of blocks in the South-North and West-East directions, respectively.
     spacing : float, tuple = (s_north, s_east), or None
         The block size in the South-North and West-East directions,
         respectively. A single value means that the size is equal in both
@@ -82,12 +88,14 @@ class BlockReduce(BaseEstimator):  # pylint: disable=too-few-public-methods
     def __init__(
         self,
         reduction,
-        spacing,
+        spacing=None,
         region=None,
         adjust="spacing",
         center_coordinates=False,
+        shape=None,
     ):
         self.reduction = reduction
+        self.shape = shape
         self.spacing = spacing
         self.region = region
         self.adjust = adjust
@@ -136,7 +144,11 @@ class BlockReduce(BaseEstimator):  # pylint: disable=too-few-public-methods
             coordinates, data, weights, unpack=False
         )
         blocks, labels = block_split(
-            coordinates, self.spacing, self.adjust, self.region
+            coordinates,
+            spacing=self.spacing,
+            shape=self.shape,
+            adjust=self.adjust,
+            region=self.region,
         )
         if any(w is None for w in weights):
             reduction = self.reduction
@@ -241,6 +253,10 @@ class BlockMean(BlockReduce):  # pylint: disable=too-few-public-methods
     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.
 
+    The size of the blocks can be specified by the *spacing* parameter. Alternatively,
+    the number of blocks in the South-North and West-East directions can be specified
+    using the *shape* parameter.
+
     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
     will be rounded to the nearest multiple. Optionally, the East and North
@@ -254,6 +270,8 @@ class BlockMean(BlockReduce):  # pylint: disable=too-few-public-methods
 
     Parameters
     ----------
+    shape : tuple = (n_north, n_east) or None
+        The number of blocks in the South-North and West-East directions, respectively.
     spacing : float, tuple = (s_north, s_east), or None
         The block size in the South-North and West-East directions,
         respectively. A single value means that the size is equal in both
@@ -287,13 +305,21 @@ class BlockMean(BlockReduce):  # pylint: disable=too-few-public-methods
 
     def __init__(
         self,
-        spacing,
+        spacing=None,
         region=None,
         adjust="spacing",
         center_coordinates=False,
         uncertainty=False,
+        shape=None,
     ):
-        super().__init__(np.average, spacing, region, adjust, center_coordinates)
+        super().__init__(
+            reduction=np.average,
+            shape=shape,
+            spacing=spacing,
+            region=region,
+            adjust=adjust,
+            center_coordinates=center_coordinates,
+        )
         self.uncertainty = uncertainty
 
     def filter(self, coordinates, data, weights=None):
@@ -342,7 +368,11 @@ class BlockMean(BlockReduce):  # pylint: disable=too-few-public-methods
                 "'uncertainty=False' to produce variance weights instead."
             )
         blocks, labels = block_split(
-            coordinates, self.spacing, self.adjust, self.region
+            coordinates,
+            spacing=self.spacing,
+            shape=self.shape,
+            adjust=self.adjust,
+            region=self.region,
         )
         ncomps = len(data)
         columns = {"data{}".format(i): comp.ravel() for i, comp in enumerate(data)}
diff --git a/verde/coordinates.py b/verde/coordinates.py
index 1903394..f64724f 100644
--- a/verde/coordinates.py
+++ b/verde/coordinates.py
@@ -658,12 +658,16 @@ def inside(coordinates, region):
     return are_inside
 
 
-def block_split(coordinates, spacing, adjust="spacing", region=None):
+def block_split(coordinates, spacing=None, adjust="spacing", region=None, shape=None):
     """
     Split a region into blocks and label points according to where they fall.
 
-    The labels are integers corresponding to the index of the block. The same
-    index is used for the coordinates of each block.
+    The labels are integers corresponding to the index of the block. Also returns the
+    coordinates of the center of each block (following the same index as the labels).
+
+    The size of the blocks can be specified by the *spacing* parameter. Alternatively,
+    the number of blocks in the South-North and West-East directions can be specified
+    using the *shape* parameter.
 
     .. note::
 
@@ -676,6 +680,8 @@ def block_split(coordinates, spacing, adjust="spacing", region=None):
         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.
+    shape : tuple = (n_north, n_east) or None
+        The number of blocks in the South-North and West-East directions, respectively.
     spacing : float, tuple = (s_north, s_east), or None
         The block size in the South-North and West-East directions,
         respectively. A single value means that the size is equal in both
@@ -719,6 +725,19 @@ def block_split(coordinates, spacing, adjust="spacing", region=None):
      [2 2 2 3 3 3]
      [2 2 2 3 3 3]
      [2 2 2 3 3 3]]
+    >>> # Use the shape instead of the block size
+    >>> block_coords, labels = block_split(coords, shape=(4, 2))
+    >>> for coord in block_coords:
+    ...     print(', '.join(['{:.3f}'.format(i) for i in coord]))
+    -3.750, -1.250, -3.750, -1.250, -3.750, -1.250, -3.750, -1.250
+    5.625, 5.625, 6.875, 6.875, 8.125, 8.125, 9.375, 9.375
+    >>> print(labels.reshape(coords[0].shape))
+    [[0 0 0 1 1 1]
+     [0 0 0 1 1 1]
+     [2 2 2 3 3 3]
+     [4 4 4 5 5 5]
+     [6 6 6 7 7 7]
+     [6 6 6 7 7 7]]
 
     """
     if region is None:
@@ -726,7 +745,7 @@ def block_split(coordinates, spacing, adjust="spacing", region=None):
     block_coords = tuple(
         i.ravel()
         for i in grid_coordinates(
-            region, spacing=spacing, adjust=adjust, pixel_register=True
+            region, spacing=spacing, shape=shape, adjust=adjust, pixel_register=True
         )
     )
     tree = kdtree(block_coords)


Test patch:
diff --git a/verde/tests/test_blockreduce.py b/verde/tests/test_blockreduce.py
index 53b8dce..66b6969 100644
--- a/verde/tests/test_blockreduce.py
+++ b/verde/tests/test_blockreduce.py
@@ -26,6 +26,21 @@ def test_block_reduce():
     npt.assert_allclose(block_coords[1][::5], np.linspace(5.5, 9.5, 5))
 
 
+def test_block_reduce_shape():
+    "Try reducing constant values in a regular grid using shape"
+    region = (-5, 0, 5, 10)
+    east, north = grid_coordinates(region, spacing=0.1, pixel_register=True)
+    data = 20 * np.ones_like(east)
+    reducer = BlockReduce(np.mean, shape=(5, 5))
+    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
+    npt.assert_allclose(block_data, 20)
+    npt.assert_allclose(block_coords[0][:5], np.linspace(-4.5, -0.5, 5))
+    npt.assert_allclose(block_coords[1][::5], np.linspace(5.5, 9.5, 5))
+
+
 def test_block_reduce_scatter():
     "Try reducing constant values in a dense enough scatter"
     region = (-5, 0, 5, 10)