From 45fdf18156f991c634c72a1a058e2264a687f3b6 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Fri, 15 May 2026 17:41:47 +0100
Subject: [PATCH 01/16] Add generate_data.py script for experiments

---
 generate_data.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100755 generate_data.py

diff --git a/generate_data.py b/generate_data.py
new file mode 100755
index 0000000..f90dfcc
--- /dev/null
+++ b/generate_data.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+
+import json
+
+import ants.io.save
+import ants.tests.stock
+import numpy as np
+
+SHAPE = (216, 432)
+
+source_data_1 = np.arange(np.prod(SHAPE), dtype=np.float64).reshape(SHAPE)
+source_cube_1 = ants.tests.stock.geodetic(data=source_data_1)
+
+source_data_2 = np.zeros_like(source_data_1)
+source_cube_2 = ants.tests.stock.geodetic(data=source_data_2)
+
+ants.io.save.netcdf(source_cube_1, ".scratch/source_1.nc", netcdf_format="NETCDF4")
+ants.io.save.netcdf(source_cube_2, ".scratch/source_2.nc", netcdf_format="NETCDF4")
+
+
+# generate polygon
+theta = np.linspace(0, np.pi * 2)
+r = 30.0 + 10 * np.sin(theta * 5)
+x = r * np.cos(theta)
+y = r * np.sin(theta)
+xy = np.stack((x, y)).T
+with open(".scratch/polygon.json", "w") as f:
+    json.dump(xy.tolist(), f)

From 2d89adb0b8a893e3e123bb2a3f7f16e62c5c1587 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Fri, 15 May 2026 17:43:11 +0100
Subject: [PATCH 02/16] Add (unused) blending_distance argument, and rename
 some variables

---
 lib/ants/analysis/__init__.py      |  4 ++--
 lib/ants/analysis/_merge.py        | 24 ++++++++++++------------
 lib/ants/cli/ancil_fill_n_merge.py |  7 ++++++-
 3 files changed, 20 insertions(+), 15 deletions(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index cddb441..3aa6148 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -187,7 +187,7 @@ def standard_deviation(source, src_mean):
     return awm
 
 
-def merge(primary_cube, alternate_cube, validity_polygon=None):
+def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance=None):
     """
     Merges data from the alternative cube into the primary cube.
 
@@ -241,7 +241,7 @@ def merge(primary_cube, alternate_cube, validity_polygon=None):
     )
     result = iris.cube.CubeList([])
     for src1, src2 in zip(primary_cubes, alternate_cubes):
-        nsource = _merge.merge(src1, src2, validity_polygon)
+        nsource = _merge.merge(src1, src2, validity_polygon, blending_distance)
         result.append(nsource)
     if isinstance(primary_cube, iris.cube.Cube):
         result = result[0]
diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 7a2d025..b61cec5 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -485,7 +485,7 @@ def _unified_grid(cube, cube2):
     return merged_cube
 
 
-def merge(primary_cube, alternate_cube, validity_polygon=None):
+def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance=None):
     """
     Merges data from the alternative cube into the primary cube.
 
@@ -630,21 +630,21 @@ def _overwrite_data(cube, cube2):
 
         # Apply data which is outside of the polygon to the other.
         # Transpose the data (view) to allow broadcasting
-        pdata, pmask = horizontal_grid_reorder(merged_cube)
-        adata, amask = horizontal_grid_reorder(full_alternate_cube)
-        pdata[full_mask_outside] = adata[full_mask_outside]
-        pmask[full_mask_outside] = amask[full_mask_outside]
+        primary_data, primary_mask = horizontal_grid_reorder(merged_cube)
+        alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
+        primary_data[full_mask_outside] = alternate_data[full_mask_outside]
+        primary_mask[full_mask_outside] = alternate_mask[full_mask_outside]
 
     # Identify overlap priority using np.nan values (these are assigned
     # when source cells are beyond the extent of the target grid whilst
     # regridding).
-    pdata, pmask = horizontal_grid_reorder(merged_cube)
-    adata, amask = horizontal_grid_reorder(full_alternate_cube)
-    slices = [slice(None)] * pdata.ndim
-    slices[2:] = [0] * (pdata.ndim - 2)
-    nan_mask = np.isnan(pdata[tuple(slices)])
-    pdata[nan_mask] = adata[nan_mask]
-    pmask[nan_mask] = amask[nan_mask]
+    primary_data, primary_mask = horizontal_grid_reorder(merged_cube)
+    alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
+    slices = [slice(None)] * primary_data.ndim
+    slices[2:] = [0] * (primary_data.ndim - 2)
+    nan_mask = np.isnan(primary_data[tuple(slices)])
+    primary_data[nan_mask] = alternate_data[nan_mask]
+    primary_mask[nan_mask] = alternate_mask[nan_mask]
 
     # Clear some memory if we can
     np.ma.MaskedArray.shrink_mask(merged_cube.data)
diff --git a/lib/ants/cli/ancil_fill_n_merge.py b/lib/ants/cli/ancil_fill_n_merge.py
index 9bb7cea..040609a 100755
--- a/lib/ants/cli/ancil_fill_n_merge.py
+++ b/lib/ants/cli/ancil_fill_n_merge.py
@@ -96,6 +96,7 @@ def main(
     end,
     netcdf_only,
     search_method,
+    blending_distance,
 ):
     """
     Perform merge and fill operation on the provided sources.
@@ -163,7 +164,9 @@ def main(
 
     result = primary_cubes
     if alternate_cubes is not None:
-        result = ants.analysis.merge(primary_cubes, alternate_cubes, validity_polygon)
+        result = ants.analysis.merge(
+            primary_cubes, alternate_cubes, validity_polygon, blending_distance
+        )
     if target_mask_filepath:
         ants.analysis.make_consistent_with_lsm(result, lbm, invert_mask, search_method)
 
@@ -226,6 +229,7 @@ def _get_parser():
         required=False,
         default="spiral",
     )
+    parser.add_argument("--blending-distance", type=float)
     return parser
 
 
@@ -251,6 +255,7 @@ def cli_interface():
         end=args.end,
         netcdf_only=args.netcdf_only,
         search_method=args.search_method,
+        blending_distance=args.blending_distance,
     )
 
 

From a2a82ef1339c8b4a93a552f84d6714941b35e393 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Fri, 15 May 2026 18:30:33 +0100
Subject: [PATCH 03/16] Add initial blending functionality

---
 lib/ants/analysis/_merge.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index b61cec5..c86dea0 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -20,6 +20,7 @@
 import numpy.lib.stride_tricks as stride
 import shapely
 from pykdtree.kdtree import KDTree
+from scipy.ndimage import distance_transform_edt
 from shapely.geometry import Polygon
 from shapely.vectorized import contains
 
@@ -632,7 +633,12 @@ def _overwrite_data(cube, cube2):
         # Transpose the data (view) to allow broadcasting
         primary_data, primary_mask = horizontal_grid_reorder(merged_cube)
         alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
-        primary_data[full_mask_outside] = alternate_data[full_mask_outside]
+        if blending_distance:
+            primary_data[...] = blend_data(
+                primary_data, alternate_data, full_mask_outside, blending_distance
+            )
+        else:
+            primary_data[full_mask_outside] = alternate_data[full_mask_outside]
         primary_mask[full_mask_outside] = alternate_mask[full_mask_outside]
 
     # Identify overlap priority using np.nan values (these are assigned
@@ -656,6 +662,13 @@ def _overwrite_data(cube, cube2):
     return merged_cube
 
 
+def blend_data(primary_data, alternate_data, full_mask_outside, blending_distance):
+    distance_outside_polygon = distance_transform_edt(full_mask_outside)
+    outside_weight = np.clip(distance_outside_polygon / blending_distance, 0.0, 1.0)
+    blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
+    return blended
+
+
 def _spiral_wrapper(
     unresolved_mask,
     to_fill_mask,

From e8cb770d5886c13ec78cfb0fd4983c0e3f226acc Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 09:39:30 +0100
Subject: [PATCH 04/16] Add documentation for blending

---
 lib/ants/analysis/__init__.py      | 17 ++++++-
 lib/ants/analysis/_merge.py        | 78 +++++++++++++++++++++++++++++-
 lib/ants/cli/ancil_fill_n_merge.py | 21 ++++++--
 3 files changed, 109 insertions(+), 7 deletions(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index 3aa6148..7254c08 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -195,8 +195,15 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
     cube which lay outside the provided polygon, override the values of the
     primary at those locations.  Containment is defined as any cell corner
     which lies within the polygon.  "Within" explicitly does not include
-    those points which exactly lay on the polygon boundary.  Where multiple
-    primary and alternate cubes are provided, then these are paired
+    those points which exactly lay on the polygon boundary.
+
+    A blending between the sources can be applied by specifying the
+    ``blending_distance`` (for no blending, pass ``None``). A linear blending
+    between the primary and alternate sources will be applied in the region
+    immediately outside the polygon over the blending distance.
+    Beyond the blending distance, the alternate source is used.
+
+    Where multiple primary and alternate cubes are provided, then these are paired
     appropriately where possible.  Where these datasets are not defined on the
     same grid, the user should consider a regrid first to then utilise merge.
 
@@ -220,6 +227,12 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         stacked together with the primary_cube taking priority over
         alternate_cube in the case of an overlap. A runtime error will be
         raised if the primary_cube is wholly within the validity_polygon.
+    blending_distance : float
+        Distance over which blending between the primary and alternate sources
+        is applied. Note that this is in units of grid cells, not a physical distance.
+        If ``None``, no blending is applied, and there will be a hard edge between
+        the two sources.
+
 
     Returns
     -------
diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index c86dea0..0fb53ec 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -500,6 +500,12 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
     in the primary cube dataset, overrides that elements validity defined in
     the cases of a specified validity polygon.
 
+    A blending between the sources can be applied by specifying the
+    ``blending_distance`` (for no blending, pass ``None``). A linear blending
+    between the primary and alternate sources will be applied in the region
+    immediately outside the polygon over the blending distance.
+    Beyond the blending distance, the alternate source is used.
+
     Parameters
     ----------
     primary_cube : `~iris.cube.Cube`
@@ -520,6 +526,11 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         alternate_cube in the case of an overlap.
         If a validity polygon is provided and the entire primary_cube dataset
         is within the polygon then a Runtime error will be raised.
+    blending_distance : float
+        Distance over which blending between the primary and alternate sources
+        is applied. Note that this is in units of grid cells, not a physical distance.
+        If ``None``, no blending is applied, and there will be a hard edge between
+        the two sources.
 
     Raises
     ------
@@ -662,8 +673,71 @@ def _overwrite_data(cube, cube2):
     return merged_cube
 
 
-def blend_data(primary_data, alternate_data, full_mask_outside, blending_distance):
-    distance_outside_polygon = distance_transform_edt(full_mask_outside)
+def blend_data(
+    primary_data: np.ndarray,
+    alternate_data: np.ndarray,
+    mask_outside: np.ndarray,
+    blending_distance: float,
+):
+    """Blend two data sources across a specified blending distance.
+
+    Returns an array with a weighted combination of data selected from the
+    primary and alternate sources, as determined by the provided mask.
+    This is calculated as follows:
+
+    1. For all points where ``mask == False``, use the primary source
+       (call this the "primary region").
+    2. For all points where ``mask == True``, determine the distance to the nearest
+       point in the primary region.
+    3. If this distance is greater than the blending distance, use the alternate source.
+    4. If this distance is less than the blending distance, weight the two datasets
+       using a linear combination: blended = w * alternate + (1 - w) * primary,
+       where w = distance / blending_distance.
+
+    The following diagram illustrates the blending in one dimension, with a
+    blending_distance of 4.
+
+    secondary              ___________
+                          /
+                         /
+    primary  ___________/
+
+    mask     0000000000011111111111111
+
+    Parameters
+    ----------
+    primary_data : np.ndarray
+        Source data to be used
+    alternate_data : np.ndarray
+        Source data to be used
+    mask_outside : np.ndarray
+        A boolean mask identifying the two regions: False for the primary source
+        region and True for the alternate source region.
+    blending_distance : float
+        Distance over which blending between the primary and alternate sources
+        is applied. Note that this is in units of grid cells, not a physical distance.
+        As such, this is resolution dependent. See notes for more detail.
+
+    Returns
+    -------
+    blended : nd.ndarray
+        The blended data
+
+    Notes
+    -----
+    The three arrays ``primary_data``, ``alternate_data`` and ``full_mask_outside``
+    must have the same shape.
+
+    This function uses :func:`scipy.ndimage.distance_transform_edt` to calculate
+    distances between points on the grid. As such, it has no knowledge of physical
+    distance or coordinate reference systems.
+
+    Warning
+    -------
+    This function does not support masked arrays. Passing a masked array may result
+    in unexpected behaviour.
+    """
+    distance_outside_polygon = distance_transform_edt(mask_outside)
     outside_weight = np.clip(distance_outside_polygon / blending_distance, 0.0, 1.0)
     blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
     return blended
diff --git a/lib/ants/cli/ancil_fill_n_merge.py b/lib/ants/cli/ancil_fill_n_merge.py
index 040609a..cf83a4d 100755
--- a/lib/ants/cli/ancil_fill_n_merge.py
+++ b/lib/ants/cli/ancil_fill_n_merge.py
@@ -105,8 +105,13 @@ def main(
     to be provided, and may optionally have a ``polygon`` shapefile.  The
     resulting data takes values from the ``primary_source`` within the ``polygon``
     (or everywhere where valid data is present, if the ``polygon`` is not
-    provided), and values from the ``alternate_source`` everywhere else.  See
-    :func:`ants.analysis.merge` for further details.
+    provided), and values from the ``alternate_source`` everywhere else.
+    A blending between the sources can be applied by specifying the
+    ``blending_distance`` (for no blending, pass ``None``). A linear blending
+    between the primary and alternate sources will be applied in the region
+    immediately outside the polygon over the blending distance.
+    Beyond the blending distance, the alternate source is used.
+    See :func:`ants.analysis.merge` for further details.
 
     The fill stage replaces missing data values with valid data, where missing
     is defined as data that is either masked or NaN.  If a landseamask is
@@ -142,6 +147,12 @@ def main(
     search_method : :obj:`str`
         Select the search method to be used when filling missing points. The methods
         currently supported are "spiral" and "kdtree".
+    blending_distance : float
+        Distance over which blending between the primary and alternate sources
+        is applied. Note that this is in units of grid cells, not a physical distance.
+        If ``None``, no blending is applied, and there will be a hard edge between
+        the two sources.
+
     Returns
     -------
     : :class:`~iris.cube.CubeList`
@@ -229,7 +240,11 @@ def _get_parser():
         required=False,
         default="spiral",
     )
-    parser.add_argument("--blending-distance", type=float)
+    blending_help = (
+        "Distance over which blending between the primary and alternate sources "
+        "is applied. Note that this is in units of grid cells, not a physical distance."
+    )
+    parser.add_argument("--blending-distance", type=float, help=blending_help)
     return parser
 
 

From 7ab6fb7f16798bc44d7a4c00f54b843feab0ec47 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 10:23:38 +0100
Subject: [PATCH 05/16] Error when blending is provided without polygon

---
 lib/ants/analysis/__init__.py         |  8 ++++++--
 lib/ants/tests/analysis/test_merge.py | 20 ++++++++++++++++++--
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index 7254c08..47fd01c 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -231,8 +231,7 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         Distance over which blending between the primary and alternate sources
         is applied. Note that this is in units of grid cells, not a physical distance.
         If ``None``, no blending is applied, and there will be a hard edge between
-        the two sources.
-
+        the two sources. This option is only valid with a provided validity polygon.
 
     Returns
     -------
@@ -245,6 +244,11 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         If the primary cube is wholly within a provided ``validity_polygon``.
 
     """
+    if blending_distance and validity_polygon is None:
+        raise ValueError(
+            "blending_distance can only be used with a validity_polygon. "
+            f"No polygon was provided, but got {blending_distance=}"
+        )
     primary_cubes = ants.utils.cube.as_cubelist(primary_cube)
     alternate_cubes = ants.utils.cube.as_cubelist(alternate_cube)
 
diff --git a/lib/ants/tests/analysis/test_merge.py b/lib/ants/tests/analysis/test_merge.py
index 1e2af77..1ad2d44 100644
--- a/lib/ants/tests/analysis/test_merge.py
+++ b/lib/ants/tests/analysis/test_merge.py
@@ -7,6 +7,7 @@
 import ants.tests
 import iris
 import numpy as np
+import pytest
 from ants.analysis import merge
 
 
@@ -40,7 +41,22 @@ def test_call_args(self):
         alternate_cube = self.generate_dummy_cube(shape=(4, 8))
 
         with mock.patch("ants.analysis._merge.merge") as mock_method:
-            merge(primary_cube, alternate_cube, None)
+            merge(primary_cube, alternate_cube, None, None)
 
-        mock_method.assert_called_once_with(primary_cube, alternate_cube, None)
+        mock_method.assert_called_once_with(primary_cube, alternate_cube, None, None)
         self.assertFalse(self.mock_fill.called)
+
+    def test_blending_distance_no_polygon(self):
+        primary_cube = self.generate_dummy_cube(shape=(4, 8))
+        alternate_cube = self.generate_dummy_cube(shape=(4, 8))
+        expected_msg = (
+            "blending_distance can only be used with a validity_polygon. "
+            "No polygon was provided, but got blending_distance=1.0"
+        )
+        with pytest.raises(ValueError, match=expected_msg):
+            merge(
+                primary_cube,
+                alternate_cube,
+                validity_polygon=None,
+                blending_distance=1.0,
+            )

From 1eb126a1c37488f76fc4919002d88d4895103b1e Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 11:06:33 +0100
Subject: [PATCH 06/16] Add tests for blend_data

---
 lib/ants/analysis/_merge.py                   |  4 ++
 .../tests/analysis/merge/test_blend_data.py   | 64 +++++++++++++++++++
 2 files changed, 68 insertions(+)
 create mode 100644 lib/ants/tests/analysis/merge/test_blend_data.py

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 0fb53ec..d2ef55e 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -737,6 +737,10 @@ def blend_data(
     This function does not support masked arrays. Passing a masked array may result
     in unexpected behaviour.
     """
+    if blending_distance <= 0:
+        raise ValueError(
+            f"Invalid blending_distance: {blending_distance}. Must be greater than zero"
+        )
     distance_outside_polygon = distance_transform_edt(mask_outside)
     outside_weight = np.clip(distance_outside_polygon / blending_distance, 0.0, 1.0)
     blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
new file mode 100644
index 0000000..d96f55c
--- /dev/null
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -0,0 +1,64 @@
+# (C) Crown Copyright, Met Office. All rights reserved.
+#
+# This file is part of ANTS and is released under the BSD 3-Clause license.
+# See LICENSE.txt in the root of the repository for full licensing details.
+import numpy as np
+import pytest
+from ants.analysis._merge import blend_data
+
+
+def test_invalid_blending_distance():
+    primary = np.array([0])
+    alternate = np.array([0])
+    mask = np.array([0])
+    blending_distance = 0
+    expected_msg = "Invalid blending_distance: 0. Must be greater than zero"
+    with pytest.raises(ValueError, match=expected_msg):
+        blend_data(primary, alternate, mask, blending_distance)
+
+
+def test_blending_1D():
+    primary = np.zeros((7,), dtype=np.float64)
+    alternate = np.ones_like(primary)
+    mask = np.array([0, 0, 1, 1, 1, 1, 1], dtype=bool)
+    blending_distance = 4
+
+    expected = np.array([0, 0, 0.25, 0.5, 0.75, 1, 1])
+
+    blended = blend_data(primary, alternate, mask, blending_distance)
+
+    np.testing.assert_array_equal(blended, expected)
+
+
+def test_blending_2D():
+    primary = np.zeros((7, 7), dtype=np.float64)
+    alternate = np.ones_like(primary)
+    mask = np.array(
+        [
+            [1, 1, 1, 1, 1, 1, 1],
+            [1, 1, 1, 1, 1, 1, 1],
+            [1, 1, 1, 1, 1, 1, 1],
+            [1, 1, 0, 0, 0, 1, 1],
+            [1, 0, 0, 0, 0, 0, 1],
+            [1, 0, 0, 0, 0, 0, 1],
+            [1, 0, 0, 0, 0, 0, 1],
+        ],
+        dtype=bool,
+    )
+    blending_distance = 2.5
+
+    expected = np.array(
+        [
+            [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
+            [1.0, 0.89442719, 0.8, 0.8, 0.8, 0.89442719, 1.0],
+            [0.89442719, 0.56568542, 0.4, 0.4, 0.4, 0.56568542, 0.89442719],
+            [0.56568542, 0.4, 0.0, 0.0, 0.0, 0.4, 0.56568542],
+            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
+            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
+            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
+        ]
+    )
+
+    blended = blend_data(primary, alternate, mask, blending_distance)
+
+    np.testing.assert_array_almost_equal(blended, expected)

From 23ca41b3cd3a5f00f5ac5efc45073bc3506ba544 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 13:36:08 +0100
Subject: [PATCH 07/16] Add is_single_level util

---
 .../tests/utils/cube/test_is_single_level.py  | 34 +++++++++++++++++++
 lib/ants/utils/cube.py                        | 21 ++++++++++++
 2 files changed, 55 insertions(+)
 create mode 100644 lib/ants/tests/utils/cube/test_is_single_level.py

diff --git a/lib/ants/tests/utils/cube/test_is_single_level.py b/lib/ants/tests/utils/cube/test_is_single_level.py
new file mode 100644
index 0000000..87515dd
--- /dev/null
+++ b/lib/ants/tests/utils/cube/test_is_single_level.py
@@ -0,0 +1,34 @@
+# (C) Crown Copyright, Met Office. All rights reserved.
+#
+# This file is part of ANTS and is released under the BSD 3-Clause license.
+# See LICENSE.txt in the root of the repository for full licensing details.
+import ants.tests.stock
+from ants.utils.cube import is_single_level
+
+
+def test_geodetic():
+    cube = ants.tests.stock.geodetic((5, 5))
+    assert is_single_level(cube) is True
+
+
+def test_geodetic_transposed():
+    cube = ants.tests.stock.geodetic((5, 5))
+    cube.transpose()
+    assert is_single_level(cube) is True
+
+
+def test_simple_4d_with_hybrid_height():
+    cube = ants.tests.stock.simple_4d_with_hybrid_height()
+    assert is_single_level(cube) is False
+
+
+def test_simple_3d_time_varying():
+    cube = ants.tests.stock.simple_3d_time_varying()
+    assert is_single_level(cube) is False
+
+
+def test_time_and_latitude():
+    # construct cube to have time and latitude coordinates only
+    cube = ants.tests.stock.simple_3d_time_varying()[..., 0]
+    assert cube.ndim == 2
+    assert is_single_level(cube) is False
diff --git a/lib/ants/utils/cube.py b/lib/ants/utils/cube.py
index b47a80e..076e919 100644
--- a/lib/ants/utils/cube.py
+++ b/lib/ants/utils/cube.py
@@ -1560,3 +1560,24 @@ def fetch_seed_index(cube, seed):
     xd = abs(x.points - seed[1]).argmin()
     yd = abs(y.points - seed[0]).argmin()
     return xd, yd
+
+
+def is_single_level(cube: iris.cube.Cube) -> bool:
+    """Determine if a cube is defined on a single horizontal level.
+
+    A cube is identified as single level if it is 2-dimensional, and those
+    dimensions correspond to the x and y axes (in any order).
+
+    Parameters
+    ----------
+    cube: iris.cube.Cube
+        The cube to check
+
+    Returns
+    -------
+    bool
+        Whether the cube is defined on a single horizontal level
+    """
+    axes = {iris.util.guess_coord_axis(coord).lower() for coord in cube.dim_coords}
+    condition = cube.ndim == 2 and axes == {"x", "y"}
+    return condition

From 2087d486f8603421e426abf238c24fd50133fadc Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 13:45:27 +0100
Subject: [PATCH 08/16] Add check for single level field when blending

---
 lib/ants/analysis/__init__.py         |  8 +++++++-
 lib/ants/tests/analysis/test_merge.py | 20 ++++++++++++++++++++
 2 files changed, 27 insertions(+), 1 deletion(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index 47fd01c..3dc7476 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -231,7 +231,8 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         Distance over which blending between the primary and alternate sources
         is applied. Note that this is in units of grid cells, not a physical distance.
         If ``None``, no blending is applied, and there will be a hard edge between
-        the two sources. This option is only valid with a provided validity polygon.
+        the two sources. This option is only valid with a provided validity polygon,
+        and with a single level field.
 
     Returns
     -------
@@ -249,6 +250,11 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
             "blending_distance can only be used with a validity_polygon. "
             f"No polygon was provided, but got {blending_distance=}"
         )
+    if blending_distance and not ants.utils.cube.is_single_level(primary_cube):
+        raise ValueError(
+            "Blending is only supported for single level data sources. "
+            "The primary data source is not single level"
+        )
     primary_cubes = ants.utils.cube.as_cubelist(primary_cube)
     alternate_cubes = ants.utils.cube.as_cubelist(alternate_cube)
 
diff --git a/lib/ants/tests/analysis/test_merge.py b/lib/ants/tests/analysis/test_merge.py
index 1ad2d44..1c6530e 100644
--- a/lib/ants/tests/analysis/test_merge.py
+++ b/lib/ants/tests/analysis/test_merge.py
@@ -60,3 +60,23 @@ def test_blending_distance_no_polygon(self):
                 validity_polygon=None,
                 blending_distance=1.0,
             )
+
+    def test_blending_distance_multi_level(self):
+        primary_cube = ants.tests.stock.simple_3d_time_varying()
+        alternate_cube = primary_cube.copy()
+
+        # doesn't matter what the validity polygon is, as long as *something* is passed
+        validity_polygon = [[0, 0], [1, 1]]
+
+        expected_msg = (
+            "Blending is only supported for single level data sources. "
+            "The primary data source is not single level"
+        )
+
+        with pytest.raises(ValueError, match=expected_msg):
+            merge(
+                primary_cube,
+                alternate_cube,
+                validity_polygon=validity_polygon,
+                blending_distance=1.0,
+            )

From 104394d129126ee3afd95f6dee890c37ae712a30 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 13:54:38 +0100
Subject: [PATCH 09/16] Add tests for different array shapes

---
 lib/ants/analysis/_merge.py                   | 12 ++++++++
 .../tests/analysis/merge/test_blend_data.py   | 28 +++++++++++++++++++
 2 files changed, 40 insertions(+)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index d2ef55e..98f872d 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -741,6 +741,18 @@ def blend_data(
         raise ValueError(
             f"Invalid blending_distance: {blending_distance}. Must be greater than zero"
         )
+    if primary_data.shape != alternate_data.shape:
+        raise ValueError(
+            "Cannot blend sources with different shapes. "
+            f"Primary shape: {primary_data.shape}, "
+            f"Alternate shape: {alternate_data.shape}"
+        )
+    if primary_data.shape != mask_outside.shape:
+        raise ValueError(
+            "Cannot blend sources as mask shape is inconsistent with source shape. "
+            f"Source shape: {primary_data.shape}, Mask shape: {mask_outside.shape}"
+        )
+
     distance_outside_polygon = distance_transform_edt(mask_outside)
     outside_weight = np.clip(distance_outside_polygon / blending_distance, 0.0, 1.0)
     blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
index d96f55c..1e036d0 100644
--- a/lib/ants/tests/analysis/merge/test_blend_data.py
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -62,3 +62,31 @@ def test_blending_2D():
     blended = blend_data(primary, alternate, mask, blending_distance)
 
     np.testing.assert_array_almost_equal(blended, expected)
+
+
+def test_different_source_shapes():
+    primary = np.zeros((2, 3))
+    alternate = np.ones((3, 2))
+    mask = np.ones_like(primary, dtype=bool)
+    blending_distance = 1
+
+    expected_msg = (
+        "Cannot blend sources with different shapes. "
+        r"Primary shape: \(2, 3\), Alternate shape: \(3, 2\)"
+    )
+    with pytest.raises(ValueError, match=expected_msg):
+        blend_data(primary, alternate, mask, blending_distance)
+
+
+def test_different_source_and_mask_shapes():
+    primary = np.zeros((2, 3))
+    alternate = np.ones_like(primary)
+    mask = np.ones((3, 2), dtype=bool)
+    blending_distance = 1
+
+    expected_msg = (
+        "Cannot blend sources as mask shape is inconsistent with source shape. "
+        r"Source shape: \(2, 3\), Mask shape: \(3, 2\)"
+    )
+    with pytest.raises(ValueError, match=expected_msg):
+        blend_data(primary, alternate, mask, blending_distance)

From d380f2b94934aedf257f0666f5bcf81086cb34ef Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Wed, 20 May 2026 17:00:55 +0100
Subject: [PATCH 10/16] More tests, and handle circular x coordinate

---
 lib/ants/analysis/_merge.py                   |  37 +++-
 .../tests/analysis/merge/test_blend_data.py   | 209 +++++++++++-------
 2 files changed, 158 insertions(+), 88 deletions(-)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 98f872d..48d7977 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -678,6 +678,7 @@ def blend_data(
     alternate_data: np.ndarray,
     mask_outside: np.ndarray,
     blending_distance: float,
+    circular: bool = False,
 ):
     """Blend two data sources across a specified blending distance.
 
@@ -737,10 +738,39 @@ def blend_data(
     This function does not support masked arrays. Passing a masked array may result
     in unexpected behaviour.
     """
+    _validate_blend_args(primary_data, alternate_data, mask_outside, blending_distance)
+
+    if circular:
+        # Pad either side of the domain to allow for wraparound in x
+        # Do not pad in y direction
+        pad_width_x = int(np.ceil(blending_distance))
+        pad_width = [(0, 0), (pad_width_x, pad_width_x)]
+        mask_outside = np.pad(mask_outside, pad_width, mode="wrap")
+
+    distance_outside_primary_region = distance_transform_edt(mask_outside)
+
+    if circular:
+        # Retrieve central slice of the padded distance array
+        npoints_x = primary_data.shape[-1]
+        slice_x = slice(pad_width_x, pad_width_x + npoints_x)
+        distance_outside_primary_region = distance_outside_primary_region[:, slice_x]
+    outside_weight = np.clip(
+        distance_outside_primary_region / blending_distance, 0.0, 1.0
+    )
+    blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
+    return blended
+
+
+def _validate_blend_args(primary_data, alternate_data, mask_outside, blending_distance):
     if blending_distance <= 0:
         raise ValueError(
             f"Invalid blending_distance: {blending_distance}. Must be greater than zero"
         )
+    if primary_data.ndim != 2:
+        raise ValueError(
+            "Can only blend 2-dimensional data, got data with "
+            f"{primary_data.ndim} dimensions"
+        )
     if primary_data.shape != alternate_data.shape:
         raise ValueError(
             "Cannot blend sources with different shapes. "
@@ -752,11 +782,8 @@ def blend_data(
             "Cannot blend sources as mask shape is inconsistent with source shape. "
             f"Source shape: {primary_data.shape}, Mask shape: {mask_outside.shape}"
         )
-
-    distance_outside_polygon = distance_transform_edt(mask_outside)
-    outside_weight = np.clip(distance_outside_polygon / blending_distance, 0.0, 1.0)
-    blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
-    return blended
+    if blending_distance > min(primary_data.shape) / 2:
+        raise ValueError("Invalid blending_distance: greater than half the domain size")
 
 
 def _spiral_wrapper(
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
index 1e036d0..616ce3e 100644
--- a/lib/ants/tests/analysis/merge/test_blend_data.py
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -7,86 +7,129 @@
 from ants.analysis._merge import blend_data
 
 
-def test_invalid_blending_distance():
-    primary = np.array([0])
-    alternate = np.array([0])
-    mask = np.array([0])
-    blending_distance = 0
-    expected_msg = "Invalid blending_distance: 0. Must be greater than zero"
-    with pytest.raises(ValueError, match=expected_msg):
-        blend_data(primary, alternate, mask, blending_distance)
-
-
-def test_blending_1D():
-    primary = np.zeros((7,), dtype=np.float64)
-    alternate = np.ones_like(primary)
-    mask = np.array([0, 0, 1, 1, 1, 1, 1], dtype=bool)
-    blending_distance = 4
-
-    expected = np.array([0, 0, 0.25, 0.5, 0.75, 1, 1])
-
-    blended = blend_data(primary, alternate, mask, blending_distance)
-
-    np.testing.assert_array_equal(blended, expected)
-
-
-def test_blending_2D():
-    primary = np.zeros((7, 7), dtype=np.float64)
-    alternate = np.ones_like(primary)
-    mask = np.array(
-        [
-            [1, 1, 1, 1, 1, 1, 1],
-            [1, 1, 1, 1, 1, 1, 1],
-            [1, 1, 1, 1, 1, 1, 1],
-            [1, 1, 0, 0, 0, 1, 1],
-            [1, 0, 0, 0, 0, 0, 1],
-            [1, 0, 0, 0, 0, 0, 1],
-            [1, 0, 0, 0, 0, 0, 1],
-        ],
-        dtype=bool,
-    )
-    blending_distance = 2.5
-
-    expected = np.array(
-        [
-            [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
-            [1.0, 0.89442719, 0.8, 0.8, 0.8, 0.89442719, 1.0],
-            [0.89442719, 0.56568542, 0.4, 0.4, 0.4, 0.56568542, 0.89442719],
-            [0.56568542, 0.4, 0.0, 0.0, 0.0, 0.4, 0.56568542],
-            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
-            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
-            [0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4],
-        ]
-    )
-
-    blended = blend_data(primary, alternate, mask, blending_distance)
-
-    np.testing.assert_array_almost_equal(blended, expected)
-
-
-def test_different_source_shapes():
-    primary = np.zeros((2, 3))
-    alternate = np.ones((3, 2))
-    mask = np.ones_like(primary, dtype=bool)
-    blending_distance = 1
-
-    expected_msg = (
-        "Cannot blend sources with different shapes. "
-        r"Primary shape: \(2, 3\), Alternate shape: \(3, 2\)"
-    )
-    with pytest.raises(ValueError, match=expected_msg):
-        blend_data(primary, alternate, mask, blending_distance)
-
-
-def test_different_source_and_mask_shapes():
-    primary = np.zeros((2, 3))
-    alternate = np.ones_like(primary)
-    mask = np.ones((3, 2), dtype=bool)
-    blending_distance = 1
-
-    expected_msg = (
-        "Cannot blend sources as mask shape is inconsistent with source shape. "
-        r"Source shape: \(2, 3\), Mask shape: \(3, 2\)"
-    )
-    with pytest.raises(ValueError, match=expected_msg):
-        blend_data(primary, alternate, mask, blending_distance)
+class TestExceptions:
+    def test_invalid_blending_distance_0(self):
+        primary = np.array([0])
+        alternate = np.array([0])
+        mask = np.array([0])
+        blending_distance = 0
+        expected_msg = "Invalid blending_distance: 0. Must be greater than zero"
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+    def test_blending_distance_too_large(self):
+        primary = np.zeros((10, 10))
+        alternate = np.ones_like(primary)
+        mask = primary.copy()
+        blending_distance = 6
+        expected_msg = "Invalid blending_distance: greater than half the domain size"
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+    def test_1D_fails(self):
+        primary = np.zeros(3)
+        alternate = np.ones_like(primary)
+        mask = primary.copy()
+        blending_distance = 2
+        expected_msg = "Can only blend 2-dimensional data, got data with 1 dimensions"
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+    def test_3D_fails(self):
+        primary = np.zeros((3, 3, 3))
+        alternate = np.ones_like(primary)
+        mask = primary.copy()
+        blending_distance = 2
+        expected_msg = "Can only blend 2-dimensional data, got data with 3 dimensions"
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+    def test_different_source_shapes(self):
+        primary = np.zeros((2, 3))
+        alternate = np.ones((3, 2))
+        mask = np.ones_like(primary, dtype=bool)
+        blending_distance = 1
+
+        expected_msg = (
+            "Cannot blend sources with different shapes. "
+            r"Primary shape: \(2, 3\), Alternate shape: \(3, 2\)"
+        )
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+    def test_different_source_and_mask_shapes(self):
+        primary = np.zeros((2, 3))
+        alternate = np.ones_like(primary)
+        mask = np.ones((3, 2), dtype=bool)
+        blending_distance = 1
+
+        expected_msg = (
+            "Cannot blend sources as mask shape is inconsistent with source shape. "
+            r"Source shape: \(2, 3\), Mask shape: \(3, 2\)"
+        )
+        with pytest.raises(ValueError, match=expected_msg):
+            blend_data(primary, alternate, mask, blending_distance)
+
+
+class TestFunctionality:
+    @pytest.fixture()
+    def primary(self):
+        return np.zeros((7, 7), dtype=np.float64)
+
+    @pytest.fixture()
+    def alternate(self):
+        return np.ones((7, 7), dtype=np.float64)
+
+    @pytest.fixture()
+    def mask(self):
+        mask = np.array(
+            [
+                [1, 1, 1, 1, 1, 1, 1],
+                [1, 1, 1, 1, 1, 1, 1],
+                [1, 1, 1, 1, 1, 1, 1],
+                [1, 1, 1, 0, 0, 1, 1],
+                [1, 1, 1, 0, 0, 0, 1],
+                [1, 1, 1, 0, 0, 0, 0],
+                [1, 1, 1, 0, 0, 0, 0],
+            ],
+            dtype=bool,
+        )
+        return mask
+
+    def test_blending_2D(self, primary, alternate, mask):
+        blending_distance = 2.5
+
+        expected = np.array(
+            [
+                [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
+                [1.0, 1.0, 0.89442719, 0.8, 0.8, 0.89442719, 1.0],
+                [1.0, 0.89442719, 0.56568542, 0.4, 0.4, 0.56568542, 0.89442719],
+                [1.0, 0.8, 0.4, 0.0, 0.0, 0.4, 0.56568542],
+                [1.0, 0.8, 0.4, 0.0, 0.0, 0.0, 0.4],
+                [1.0, 0.8, 0.4, 0.0, 0.0, 0.0, 0.0],
+                [1.0, 0.8, 0.4, 0.0, 0.0, 0.0, 0.0],
+            ]
+        )
+
+        blended = blend_data(primary, alternate, mask, blending_distance)
+
+        np.testing.assert_array_almost_equal(blended, expected)
+
+    def test_blending_2D_circular(self, primary, alternate, mask):
+        blending_distance = 2.5
+
+        expected = np.array(
+            [
+                [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
+                [1.0, 1.0, 0.89442719, 0.8, 0.8, 0.89442719, 1.0],
+                [1.0, 0.89442719, 0.56568542, 0.4, 0.4, 0.56568542, 0.89442719],
+                [0.89442719, 0.8, 0.4, 0.0, 0.0, 0.4, 0.56568542],
+                [0.56568542, 0.8, 0.4, 0.0, 0.0, 0.0, 0.4],
+                [0.4, 0.8, 0.4, 0.0, 0.0, 0.0, 0.0],
+                [0.4, 0.8, 0.4, 0.0, 0.0, 0.0, 0.0],
+            ]
+        )
+
+        blended = blend_data(primary, alternate, mask, blending_distance, circular=True)
+
+        np.testing.assert_array_almost_equal(blended, expected)

From 4a2fee3f4956c944f68fb1c8b36c065c9028bc2d Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 09:28:16 +0100
Subject: [PATCH 11/16] Pass circular flag to blend_data, refactor validation

---
 lib/ants/analysis/__init__.py | 42 ++++++++++++++++++++++++++---------
 lib/ants/analysis/_merge.py   |  7 +++++-
 2 files changed, 38 insertions(+), 11 deletions(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index 3dc7476..30577c1 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -245,19 +245,14 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
         If the primary cube is wholly within a provided ``validity_polygon``.
 
     """
-    if blending_distance and validity_polygon is None:
-        raise ValueError(
-            "blending_distance can only be used with a validity_polygon. "
-            f"No polygon was provided, but got {blending_distance=}"
-        )
-    if blending_distance and not ants.utils.cube.is_single_level(primary_cube):
-        raise ValueError(
-            "Blending is only supported for single level data sources. "
-            "The primary data source is not single level"
-        )
     primary_cubes = ants.utils.cube.as_cubelist(primary_cube)
     alternate_cubes = ants.utils.cube.as_cubelist(alternate_cube)
 
+    if blending_distance:
+        _validate_args_with_blending(
+            primary_cubes, alternate_cubes, validity_polygon, blending_distance
+        )
+
     # Group (sort) cubes so they are ordered in a way suitable for merging.
     primary_cubes, alternate_cubes = ants.utils.cube.sort_cubes(
         primary_cubes, alternate_cubes
@@ -271,6 +266,33 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
     return result
 
 
+def _validate_args_with_blending(
+    primary_cubes, alternate_cubes, validity_polygon, blending_distance
+):
+    """Specific validation for merge arguments when blending is provided."""
+    if validity_polygon is None:
+        raise ValueError(
+            "blending_distance can only be used with a validity_polygon. "
+            f"No polygon was provided, but got {blending_distance=}"
+        )
+
+    all_primary_single_level = all(map(ants.utils.cube.is_single_level, primary_cubes))
+    if not all_primary_single_level:
+        raise ValueError(
+            "Blending is only supported for single level data sources. "
+            "The primary data source is not single level"
+        )
+
+    all_alternate_single_level = all(
+        map(ants.utils.cube.is_single_level, alternate_cubes)
+    )
+    if not all_alternate_single_level:
+        raise ValueError(
+            "Blending is only supported for single level data sources. "
+            "The alternate data source is not single level"
+        )
+
+
 def _flood_fill_neighbour_identify(
     shape, coords, seed_point, extended_neighbourhood, wraparound
 ):
diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 48d7977..17f9037 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -645,8 +645,13 @@ def _overwrite_data(cube, cube2):
         primary_data, primary_mask = horizontal_grid_reorder(merged_cube)
         alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
         if blending_distance:
+            is_circular = primary_cube.coord(axis="x").circular
             primary_data[...] = blend_data(
-                primary_data, alternate_data, full_mask_outside, blending_distance
+                primary_data,
+                alternate_data,
+                full_mask_outside,
+                blending_distance,
+                is_circular,
             )
         else:
             primary_data[full_mask_outside] = alternate_data[full_mask_outside]

From 0c75ecfc9b2753f93cbc120c213b0d7b25faaf8d Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 11:30:57 +0100
Subject: [PATCH 12/16] Rename blend args for clarity

---
 lib/ants/analysis/_merge.py                   | 76 +++++++++----------
 .../tests/analysis/merge/test_blend_data.py   |  3 +-
 2 files changed, 39 insertions(+), 40 deletions(-)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 17f9037..e61ed5e 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -646,6 +646,7 @@ def _overwrite_data(cube, cube2):
         alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
         if blending_distance:
             is_circular = primary_cube.coord(axis="x").circular
+            breakpoint()
             primary_data[...] = blend_data(
                 primary_data,
                 alternate_data,
@@ -679,46 +680,46 @@ def _overwrite_data(cube, cube2):
 
 
 def blend_data(
-    primary_data: np.ndarray,
-    alternate_data: np.ndarray,
-    mask_outside: np.ndarray,
+    from_array: np.ndarray,
+    into_array: np.ndarray,
+    mask: np.ndarray,
     blending_distance: float,
     circular: bool = False,
 ):
     """Blend two data sources across a specified blending distance.
 
-    Returns an array with a weighted combination of data selected from the
-    primary and alternate sources, as determined by the provided mask.
+    Returns an array with a weighted combination of data selected from the two
+    sources, as determined by the provided mask.
+
     This is calculated as follows:
 
-    1. For all points where ``mask == False``, use the primary source
-       (call this the "primary region").
+    1. For all points where ``mask == False``, use the "from" source
     2. For all points where ``mask == True``, determine the distance to the nearest
-       point in the primary region.
-    3. If this distance is greater than the blending distance, use the alternate source.
+       point in the "from" region.
+    3. If this distance is greater than the blending distance, use the "into" source.
     4. If this distance is less than the blending distance, weight the two datasets
-       using a linear combination: blended = w * alternate + (1 - w) * primary,
+       using a linear combination: blended = w * from_array + (1 - w) * into_array,
        where w = distance / blending_distance.
 
     The following diagram illustrates the blending in one dimension, with a
     blending_distance of 4.
 
-    secondary              ___________
-                          /
-                         /
-    primary  ___________/
+    into                ___________
+                       /
+                      /
+    from  ___________/
 
-    mask     0000000000011111111111111
+    mask  0000000000011111111111111
 
     Parameters
     ----------
-    primary_data : np.ndarray
-        Source data to be used
-    alternate_data : np.ndarray
-        Source data to be used
-    mask_outside : np.ndarray
-        A boolean mask identifying the two regions: False for the primary source
-        region and True for the alternate source region.
+    from_array : np.ndarray
+        Source data to be blended from
+    into_array : np.ndarray
+        Source data to be blended into
+    mask : np.ndarray
+        A boolean mask identifying the two regions: False for the "from" source
+        region and True for the "into" source region.
     blending_distance : float
         Distance over which blending between the primary and alternate sources
         is applied. Note that this is in units of grid cells, not a physical distance.
@@ -731,7 +732,7 @@ def blend_data(
 
     Notes
     -----
-    The three arrays ``primary_data``, ``alternate_data`` and ``full_mask_outside``
+    The three arrays ``from_array``, ``into_array`` and ``mask``
     must have the same shape.
 
     This function uses :func:`scipy.ndimage.distance_transform_edt` to calculate
@@ -743,51 +744,50 @@ def blend_data(
     This function does not support masked arrays. Passing a masked array may result
     in unexpected behaviour.
     """
-    _validate_blend_args(primary_data, alternate_data, mask_outside, blending_distance)
+    _validate_blend_args(from_array, into_array, mask, blending_distance)
 
     if circular:
         # Pad either side of the domain to allow for wraparound in x
         # Do not pad in y direction
         pad_width_x = int(np.ceil(blending_distance))
         pad_width = [(0, 0), (pad_width_x, pad_width_x)]
-        mask_outside = np.pad(mask_outside, pad_width, mode="wrap")
+        mask = np.pad(mask, pad_width, mode="wrap")
 
-    distance_outside_primary_region = distance_transform_edt(mask_outside)
+    distance_outside_primary_region = distance_transform_edt(mask)
 
     if circular:
         # Retrieve central slice of the padded distance array
-        npoints_x = primary_data.shape[-1]
+        npoints_x = from_array.shape[-1]
         slice_x = slice(pad_width_x, pad_width_x + npoints_x)
         distance_outside_primary_region = distance_outside_primary_region[:, slice_x]
     outside_weight = np.clip(
         distance_outside_primary_region / blending_distance, 0.0, 1.0
     )
-    blended = (outside_weight * alternate_data) + (1 - outside_weight) * primary_data
+    blended = (outside_weight * into_array) + (1 - outside_weight) * from_array
     return blended
 
 
-def _validate_blend_args(primary_data, alternate_data, mask_outside, blending_distance):
+def _validate_blend_args(from_array, into_array, mask, blending_distance):
     if blending_distance <= 0:
         raise ValueError(
             f"Invalid blending_distance: {blending_distance}. Must be greater than zero"
         )
-    if primary_data.ndim != 2:
+    if from_array.ndim != 2:
         raise ValueError(
             "Can only blend 2-dimensional data, got data with "
-            f"{primary_data.ndim} dimensions"
+            f"{from_array.ndim} dimensions"
         )
-    if primary_data.shape != alternate_data.shape:
+    if from_array.shape != into_array.shape:
         raise ValueError(
-            "Cannot blend sources with different shapes. "
-            f"Primary shape: {primary_data.shape}, "
-            f"Alternate shape: {alternate_data.shape}"
+            f"Cannot blend sources with different shapes: "
+            f"{from_array.shape} and {into_array.shape}"
         )
-    if primary_data.shape != mask_outside.shape:
+    if from_array.shape != mask.shape:
         raise ValueError(
             "Cannot blend sources as mask shape is inconsistent with source shape. "
-            f"Source shape: {primary_data.shape}, Mask shape: {mask_outside.shape}"
+            f"Source shape: {from_array.shape}, Mask shape: {mask.shape}"
         )
-    if blending_distance > min(primary_data.shape) / 2:
+    if blending_distance > min(from_array.shape) / 2:
         raise ValueError("Invalid blending_distance: greater than half the domain size")
 
 
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
index 616ce3e..6b7820a 100644
--- a/lib/ants/tests/analysis/merge/test_blend_data.py
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -51,8 +51,7 @@ def test_different_source_shapes(self):
         blending_distance = 1
 
         expected_msg = (
-            "Cannot blend sources with different shapes. "
-            r"Primary shape: \(2, 3\), Alternate shape: \(3, 2\)"
+            r"Cannot blend sources with different shapes: \(2, 3\) and \(3, 2\)"
         )
         with pytest.raises(ValueError, match=expected_msg):
             blend_data(primary, alternate, mask, blending_distance)

From 8a34bbf154c154e194484efde97e344c0f186c8a Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 13:53:12 +0100
Subject: [PATCH 13/16] Further renaming and add warning if region is fully
 blended

---
 lib/ants/analysis/_merge.py                   | 36 ++++++-----
 .../tests/analysis/merge/test_blend_data.py   | 61 ++++++++++---------
 2 files changed, 54 insertions(+), 43 deletions(-)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index e61ed5e..fbfd3ac 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -646,13 +646,12 @@ def _overwrite_data(cube, cube2):
         alternate_data, alternate_mask = horizontal_grid_reorder(full_alternate_cube)
         if blending_distance:
             is_circular = primary_cube.coord(axis="x").circular
-            breakpoint()
             primary_data[...] = blend_data(
-                primary_data,
-                alternate_data,
-                full_mask_outside,
-                blending_distance,
-                is_circular,
+                from_array=alternate_data,
+                into_array=primary_data,
+                mask=~full_mask_outside,
+                blending_distance=blending_distance,
+                circular=is_circular,
             )
         else:
             primary_data[full_mask_outside] = alternate_data[full_mask_outside]
@@ -721,7 +720,7 @@ def blend_data(
         A boolean mask identifying the two regions: False for the "from" source
         region and True for the "into" source region.
     blending_distance : float
-        Distance over which blending between the primary and alternate sources
+        Distance over which blending between the sources
         is applied. Note that this is in units of grid cells, not a physical distance.
         As such, this is resolution dependent. See notes for more detail.
 
@@ -753,17 +752,23 @@ def blend_data(
         pad_width = [(0, 0), (pad_width_x, pad_width_x)]
         mask = np.pad(mask, pad_width, mode="wrap")
 
-    distance_outside_primary_region = distance_transform_edt(mask)
+    distance_into_region = distance_transform_edt(mask)
+    max_distance_into_region = distance_into_region.max()
+    if max_distance_into_region < blending_distance:
+        warnings.warn(
+            "All points within the blending region are within the blending distance. "
+            f"Specified {blending_distance=}, maximum distance into domain: "
+            f"{max_distance_into_region}"
+        )
 
     if circular:
         # Retrieve central slice of the padded distance array
         npoints_x = from_array.shape[-1]
         slice_x = slice(pad_width_x, pad_width_x + npoints_x)
-        distance_outside_primary_region = distance_outside_primary_region[:, slice_x]
-    outside_weight = np.clip(
-        distance_outside_primary_region / blending_distance, 0.0, 1.0
-    )
-    blended = (outside_weight * into_array) + (1 - outside_weight) * from_array
+        distance_into_region = distance_into_region[:, slice_x]
+
+    into_weight = np.clip(distance_into_region / blending_distance, 0.0, 1.0)
+    blended = (into_weight * into_array) + (1 - into_weight) * from_array
     return blended
 
 
@@ -788,7 +793,10 @@ def _validate_blend_args(from_array, into_array, mask, blending_distance):
             f"Source shape: {from_array.shape}, Mask shape: {mask.shape}"
         )
     if blending_distance > min(from_array.shape) / 2:
-        raise ValueError("Invalid blending_distance: greater than half the domain size")
+        raise ValueError(
+            f"Invalid {blending_distance=}: greater than half the domain size "
+            f"(shape={from_array.shape})"
+        )
 
 
 def _spiral_wrapper(
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
index 6b7820a..922d5fa 100644
--- a/lib/ants/tests/analysis/merge/test_blend_data.py
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -9,56 +9,59 @@
 
 class TestExceptions:
     def test_invalid_blending_distance_0(self):
-        primary = np.array([0])
-        alternate = np.array([0])
+        source1 = np.array([0])
+        source2 = np.array([0])
         mask = np.array([0])
         blending_distance = 0
         expected_msg = "Invalid blending_distance: 0. Must be greater than zero"
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
     def test_blending_distance_too_large(self):
-        primary = np.zeros((10, 10))
-        alternate = np.ones_like(primary)
-        mask = primary.copy()
+        source1 = np.zeros((10, 10))
+        source2 = np.ones_like(source1)
+        mask = source1.copy()
         blending_distance = 6
-        expected_msg = "Invalid blending_distance: greater than half the domain size"
+        expected_msg = (
+            "Invalid blending_distance=6: greater than half the domain size "
+            r"\(shape=\(10, 10\)\)"
+        )
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
     def test_1D_fails(self):
-        primary = np.zeros(3)
-        alternate = np.ones_like(primary)
-        mask = primary.copy()
+        source1 = np.zeros(3)
+        source2 = np.ones_like(source1)
+        mask = source1.copy()
         blending_distance = 2
         expected_msg = "Can only blend 2-dimensional data, got data with 1 dimensions"
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
     def test_3D_fails(self):
-        primary = np.zeros((3, 3, 3))
-        alternate = np.ones_like(primary)
-        mask = primary.copy()
+        source1 = np.zeros((3, 3, 3))
+        source2 = np.ones_like(source1)
+        mask = source1.copy()
         blending_distance = 2
         expected_msg = "Can only blend 2-dimensional data, got data with 3 dimensions"
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
     def test_different_source_shapes(self):
-        primary = np.zeros((2, 3))
-        alternate = np.ones((3, 2))
-        mask = np.ones_like(primary, dtype=bool)
+        source1 = np.zeros((2, 3))
+        source2 = np.ones((3, 2))
+        mask = np.ones_like(source1, dtype=bool)
         blending_distance = 1
 
         expected_msg = (
             r"Cannot blend sources with different shapes: \(2, 3\) and \(3, 2\)"
         )
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
     def test_different_source_and_mask_shapes(self):
-        primary = np.zeros((2, 3))
-        alternate = np.ones_like(primary)
+        source1 = np.zeros((2, 3))
+        source2 = np.ones_like(source1)
         mask = np.ones((3, 2), dtype=bool)
         blending_distance = 1
 
@@ -67,16 +70,16 @@ def test_different_source_and_mask_shapes(self):
             r"Source shape: \(2, 3\), Mask shape: \(3, 2\)"
         )
         with pytest.raises(ValueError, match=expected_msg):
-            blend_data(primary, alternate, mask, blending_distance)
+            blend_data(source1, source2, mask, blending_distance)
 
 
 class TestFunctionality:
     @pytest.fixture()
-    def primary(self):
+    def source1(self):
         return np.zeros((7, 7), dtype=np.float64)
 
     @pytest.fixture()
-    def alternate(self):
+    def source2(self):
         return np.ones((7, 7), dtype=np.float64)
 
     @pytest.fixture()
@@ -95,7 +98,7 @@ def mask(self):
         )
         return mask
 
-    def test_blending_2D(self, primary, alternate, mask):
+    def test_blending_2D(self, source1, source2, mask):
         blending_distance = 2.5
 
         expected = np.array(
@@ -110,11 +113,11 @@ def test_blending_2D(self, primary, alternate, mask):
             ]
         )
 
-        blended = blend_data(primary, alternate, mask, blending_distance)
+        blended = blend_data(source1, source2, mask, blending_distance)
 
         np.testing.assert_array_almost_equal(blended, expected)
 
-    def test_blending_2D_circular(self, primary, alternate, mask):
+    def test_blending_2D_circular(self, source1, source2, mask):
         blending_distance = 2.5
 
         expected = np.array(
@@ -129,6 +132,6 @@ def test_blending_2D_circular(self, primary, alternate, mask):
             ]
         )
 
-        blended = blend_data(primary, alternate, mask, blending_distance, circular=True)
+        blended = blend_data(source1, source2, mask, blending_distance, circular=True)
 
         np.testing.assert_array_almost_equal(blended, expected)

From 1da6caa19ff7ed62e1262949fe0d1f6ccfb852b9 Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 14:12:28 +0100
Subject: [PATCH 14/16] Add test for warning

---
 lib/ants/analysis/_merge.py                   |  2 +-
 .../tests/analysis/merge/test_blend_data.py   | 21 +++++++++++++++++++
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index fbfd3ac..39a26db 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -756,7 +756,7 @@ def blend_data(
     max_distance_into_region = distance_into_region.max()
     if max_distance_into_region < blending_distance:
         warnings.warn(
-            "All points within the blending region are within the blending distance. "
+            "All points within the region are within the blending distance. "
             f"Specified {blending_distance=}, maximum distance into domain: "
             f"{max_distance_into_region}"
         )
diff --git a/lib/ants/tests/analysis/merge/test_blend_data.py b/lib/ants/tests/analysis/merge/test_blend_data.py
index 922d5fa..2606b1e 100644
--- a/lib/ants/tests/analysis/merge/test_blend_data.py
+++ b/lib/ants/tests/analysis/merge/test_blend_data.py
@@ -72,6 +72,27 @@ def test_different_source_and_mask_shapes(self):
         with pytest.raises(ValueError, match=expected_msg):
             blend_data(source1, source2, mask, blending_distance)
 
+    def test_blending_covers_entire_region(self):
+        source1 = np.zeros((5, 5))
+        source2 = np.ones_like(source1)
+        mask = np.array(
+            [
+                [0, 0, 0, 0, 0],
+                [0, 1, 1, 1, 0],
+                [0, 1, 1, 1, 0],
+                [0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0],
+            ],
+            dtype=bool,
+        )
+        blending_distance = 2.0
+        expected_msg = (
+            "All points within the region are within the blending distance. "
+            "Specified blending_distance=2.0, maximum distance into domain: 1.0"
+        )
+        with pytest.warns(match=expected_msg):
+            blend_data(source1, source2, mask, blending_distance)
+
 
 class TestFunctionality:
     @pytest.fixture()

From 672b914f515cab587f1973624632d5051d109d7f Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 14:14:48 +0100
Subject: [PATCH 15/16] Remove synthetic test data script

---
 generate_data.py | 28 ----------------------------
 1 file changed, 28 deletions(-)
 delete mode 100755 generate_data.py

diff --git a/generate_data.py b/generate_data.py
deleted file mode 100755
index f90dfcc..0000000
--- a/generate_data.py
+++ /dev/null
@@ -1,28 +0,0 @@
-#!/usr/bin/env python
-
-import json
-
-import ants.io.save
-import ants.tests.stock
-import numpy as np
-
-SHAPE = (216, 432)
-
-source_data_1 = np.arange(np.prod(SHAPE), dtype=np.float64).reshape(SHAPE)
-source_cube_1 = ants.tests.stock.geodetic(data=source_data_1)
-
-source_data_2 = np.zeros_like(source_data_1)
-source_cube_2 = ants.tests.stock.geodetic(data=source_data_2)
-
-ants.io.save.netcdf(source_cube_1, ".scratch/source_1.nc", netcdf_format="NETCDF4")
-ants.io.save.netcdf(source_cube_2, ".scratch/source_2.nc", netcdf_format="NETCDF4")
-
-
-# generate polygon
-theta = np.linspace(0, np.pi * 2)
-r = 30.0 + 10 * np.sin(theta * 5)
-x = r * np.cos(theta)
-y = r * np.sin(theta)
-xy = np.stack((x, y)).T
-with open(".scratch/polygon.json", "w") as f:
-    json.dump(xy.tolist(), f)

From 061524ef35b4573385f5a3250418ec34e010d27c Mon Sep 17 00:00:00 2001
From: Josh Rackham <144251043+jrackham-mo@users.noreply.github.com>
Date: Thu, 21 May 2026 14:27:08 +0100
Subject: [PATCH 16/16] Update wording for new blend direction

---
 lib/ants/analysis/__init__.py | 2 +-
 lib/ants/analysis/_merge.py   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/ants/analysis/__init__.py b/lib/ants/analysis/__init__.py
index 30577c1..4000a94 100644
--- a/lib/ants/analysis/__init__.py
+++ b/lib/ants/analysis/__init__.py
@@ -200,7 +200,7 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
     A blending between the sources can be applied by specifying the
     ``blending_distance`` (for no blending, pass ``None``). A linear blending
     between the primary and alternate sources will be applied in the region
-    immediately outside the polygon over the blending distance.
+    immediately inside the polygon over the blending distance.
     Beyond the blending distance, the alternate source is used.
 
     Where multiple primary and alternate cubes are provided, then these are paired
diff --git a/lib/ants/analysis/_merge.py b/lib/ants/analysis/_merge.py
index 39a26db..704f385 100644
--- a/lib/ants/analysis/_merge.py
+++ b/lib/ants/analysis/_merge.py
@@ -503,7 +503,7 @@ def merge(primary_cube, alternate_cube, validity_polygon=None, blending_distance
     A blending between the sources can be applied by specifying the
     ``blending_distance`` (for no blending, pass ``None``). A linear blending
     between the primary and alternate sources will be applied in the region
-    immediately outside the polygon over the blending distance.
+    immediately inside the polygon over the blending distance.
     Beyond the blending distance, the alternate source is used.
 
     Parameters