From 08239e2c96133235ff37ed91833f5926e2b145dc Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Fri, 27 Oct 2017 16:11:51 +0100
Subject: [PATCH 01/10] ENH: Added conservative interpolation approach

---
 setup.py                               |   3 +
 stratify/__init__.py                   |   2 +
 stratify/_bounded_vinterp.py           | 133 ++++++++++++++++++
 stratify/_conservative.pyx             |  55 ++++++++
 stratify/tests/test_bounded_vinterp.py | 183 +++++++++++++++++++++++++
 5 files changed, 376 insertions(+)
 create mode 100644 stratify/_bounded_vinterp.py
 create mode 100644 stratify/_conservative.pyx
 create mode 100644 stratify/tests/test_bounded_vinterp.py

diff --git a/setup.py b/setup.py
index 91546ba..1ffe5ba 100644
--- a/setup.py
+++ b/setup.py
@@ -13,6 +13,9 @@
 
 extensions = [Extension('{}._vinterp'.format(NAME),
                         [os.path.join(NAME, '_vinterp.pyx')],
+                        include_dirs=[np.get_include()]),
+              Extension('{}._conservative'.format(NAME),
+                        [os.path.join(NAME, '_conservative.pyx')],
                         include_dirs=[np.get_include()])]
 
 
diff --git a/stratify/__init__.py b/stratify/__init__.py
index 9e298e6..c1af8c4 100644
--- a/stratify/__init__.py
+++ b/stratify/__init__.py
@@ -5,5 +5,7 @@
                        EXTRAPOLATE_NAN, EXTRAPOLATE_NEAREST,
                        EXTRAPOLATE_LINEAR, PyFuncExtrapolator, 
                        PyFuncInterpolator)
+from ._bounded_vinterp import interpolate as interpolate_conservative
+
 
 __version__ = '0.1a3.dev0'
diff --git a/stratify/_bounded_vinterp.py b/stratify/_bounded_vinterp.py
new file mode 100644
index 0000000..74a5e14
--- /dev/null
+++ b/stratify/_bounded_vinterp.py
@@ -0,0 +1,133 @@
+from __future__ import absolute_import
+
+import numpy as np
+
+from ._conservative import conservative_interpolation
+
+
+def interpolate(z_target, z_src, fz_src, axis=-1):
+    """
+    Interface for optimised 1d interpolation across multiple dimensions.
+
+    This function provides the ability to perform 1d interpolation on datasets
+    with more than one dimension. For instance, this function can be used to
+    interpolate a set of vertical levels, even if the interpolation coordinate
+    depends upon other dimensions.
+
+    A good use case might be when wanting to interpolate at a specific height
+    for height data which also depends on x and y - e.g. extract 1000hPa level
+    from a 3d dataset and associated pressure field. In the case of this
+    example, pressure would be the `z` coordinate, and the dataset
+    (e.g. geopotential height / temperature etc.) would be `f(z)`.
+
+    Parameters
+    ----------
+    z_target: :class:`np.ndarray`
+        Target coordinate.
+        This coordinate defines the levels to interpolate the source data
+        ``fz_src`` to.  ``z_target`` must have the same dimensionality as the
+        source coordinate ``z_src``, and the shape of ``z_target`` must match
+        the shape of ``z_src``, although the axis of interpolation may differ
+        in dimension size.
+    z_src: :class:`np.ndarray`
+        Source coordinate.
+        This coordinate defines the levels that the source data ``fz_src`` is
+        interpolated from.
+    fz_src: :class:`np.ndarray`
+        The source data; the phenomenon data values to be interpolated from
+        ``z_src`` to ``z_target``.
+        The data array must be at least ``z_src.ndim``, and its trailing
+        dimensions (i.e. those on its right hand side) must be exactly
+        the same as the shape of ``z_src``.
+    axis: int (default -1)
+        The ``fz_src`` axis to perform the interpolation over.
+
+    Returns
+    -------
+    : :class:`np.ndarray`
+        fz_src interpolated from z_src to z_target.
+
+    Note
+    ----
+    Support for 1D z_target and corresponding ND z_src will be provided in
+    future as driven by user requirement.
+
+    """
+    if z_src.ndim != z_target.ndim:
+        msg = ('Expecting source and target levels dimensionality to be '
+               'identical.  {} != {}.')
+        raise ValueError(msg.format(z_src.ndim, z_target.ndim))
+
+    # Relative axis
+    axis = axis % fz_src.ndim
+    axis_relative = axis - (fz_src.ndim - (z_target.ndim-1))
+
+    src_shape = list(z_src.shape)
+    src_shape.pop(axis_relative)
+    tgt_shape = list(z_target.shape)
+    tgt_shape.pop(axis_relative)
+
+    if src_shape != tgt_shape:
+        src_shape = list(z_src.shape)
+        src_shape[axis_relative] = '-'
+        tgt_shape = list(z_target.shape)
+        src_shape[axis_relative] = '-'
+        msg = ('Expecting the shape of the source and target levels except '
+               'the axis of interpolation to be identical.  {} != {}')
+        raise ValueError(msg.format(tuple(src_shape), tuple(tgt_shape)))
+
+    dat_shape = list(fz_src.shape)
+    dat_shape = dat_shape[-(z_src.ndim-1):]
+    src_shape = list(z_src.shape[:-1])
+    if dat_shape != src_shape:
+        dat_shape = list(fz_src.shape)
+        dat_shape[:-(z_src.ndim-1)] = '-'
+        msg = ('The provided data is not of compatible shape with the '
+               'provided source bounds. {} != {}')
+        raise ValueError(msg.format(tuple(dat_shape), tuple(src_shape)))
+
+    if z_src.shape[-1] != 2:
+        msg = 'Unexpected source and target bounds shape. shape[-1] != 2'
+        raise ValueError(msg)
+
+    # Define our source in a consistent way.
+    # [broadcasting_dims, axis_interpolation, z_varying]
+
+    # src_data
+    bdims = range(fz_src.ndim - (z_src.ndim-1))
+    data_vdims = [ind for ind in range(fz_src.ndim) if ind not in
+                  (bdims + [axis])]
+    data_transpose = bdims + [axis] + data_vdims
+    fz_src_reshaped = np.transpose(fz_src, data_transpose)
+    fz_src_orig = list(fz_src_reshaped.shape)
+    shape = (
+        int(np.product(fz_src_reshaped.shape[:len(bdims)])),
+        fz_src_reshaped.shape[len(bdims)],
+        int(np.product(fz_src_reshaped.shape[len(bdims)+1:])))
+    fz_src_reshaped = fz_src_reshaped.reshape(shape)
+
+    # Define our src and target bounds in a consistent way.
+    # [axis_interpolation, z_varying, 2]
+    vdims = list(set(range(z_src.ndim)) - set([axis_relative]))
+    z_src_reshaped = np.transpose(z_src, [axis_relative] + vdims)
+    z_target_reshaped = np.transpose(z_target, [axis_relative] + vdims)
+
+    shape = int(np.product(z_src_reshaped.shape[1:-1]))
+    z_src_reshaped = z_src_reshaped.reshape([z_src_reshaped.shape[0], shape,
+                                             z_src_reshaped.shape[-1]])
+    shape = int(np.product(z_target_reshaped.shape[1:-1]))
+    z_target_reshaped = z_target_reshaped.reshape(
+            [z_target_reshaped.shape[0], shape, z_target_reshaped.shape[-1]])
+
+    result = conservative_interpolation(
+        z_src_reshaped, z_target_reshaped, fz_src_reshaped)
+
+    # Turn the result into a shape consistent with the source.
+    # First reshape, then reverse transpose.
+    shape = fz_src_orig
+    shape[len(bdims)] = z_target.shape[axis_relative]
+    result = result.reshape(shape)
+    invert_transpose = [data_transpose.index(ind) for ind in
+                        list(range(result.ndim))]
+    result = result.transpose(invert_transpose)
+    return result
diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
new file mode 100644
index 0000000..7d8069b
--- /dev/null
+++ b/stratify/_conservative.pyx
@@ -0,0 +1,55 @@
+import numpy as np
+cimport numpy as np
+
+
+cdef calculate_weights(np.ndarray[np.float64_t, ndim=2] src_point,
+                       np.ndarray[np.float64_t, ndim=2] tgt_point):
+    # src_min    src_max
+    #    |----------|                 : Source
+    #       tgt_min     tgt_max
+    #          |------------|         : Target
+    #    |----------|                 : Delta (src_max - src_min)
+    #          |----|                 : Overlap (between src & tgt)
+    #    weight = overlap / delta
+    cdef Py_ssize_t src_ind, tgt_ind
+    cdef np.float64_t delta, weight
+    cdef np.ndarray[np.float64_t, ndim = 2] weights
+    cdef np.ndarray[np.float64_t, ndim = 1] src_cell, tgt_cell
+
+    weights = np.zeros([src_point.shape[0], tgt_point.shape[0]],
+                       dtype=np.float64)
+    for src_ind, src_cell in enumerate(src_point):
+        delta = src_cell.max() - src_cell.min()
+        for tgt_ind, tgt_cell in enumerate(tgt_point):
+            weight = (min(src_cell.max(), tgt_cell.max()) -
+                      max(src_cell.min(), tgt_cell.min())) / float(delta)
+            if weight > 0:
+                weights[src_ind, tgt_ind] = weight
+    return weights
+
+
+cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
+                   np.ndarray[np.float64_t, ndim=3] tgt_point,
+                   np.ndarray[np.float64_t, ndim=3] src_data):
+    cdef Py_ssize_t ind
+    cdef np.ndarray[np.float64_t, ndim = 3] results
+    cdef np.ndarray[np.float64_t, ndim = 2] weights
+    results = np.zeros(
+        [src_data.shape[0], tgt_point.shape[0], src_data.shape[2]],
+        dtype='float64')
+    # Calculate and apply weights
+    for ind in range(src_data.shape[2]):
+        weights = calculate_weights(src_point[:, ind], tgt_point[:, ind])
+        if not (weights.sum(axis=1) == 1).all():
+            msg = ('Weights calculation yields a less than conservative '
+                   'result.  Aborting.')
+            raise ValueError(msg)
+        results[..., ind] = (
+            weights * src_data[..., ind][..., None]).sum(axis=1)
+    return results
+
+
+def conservative_interpolation(src_points, tgt_points, src_data):
+    return apply_weights(src_points.astype('float64'),
+                         tgt_points.astype('float64'),
+                         src_data.astype('float64'))
diff --git a/stratify/tests/test_bounded_vinterp.py b/stratify/tests/test_bounded_vinterp.py
new file mode 100644
index 0000000..b7d61e9
--- /dev/null
+++ b/stratify/tests/test_bounded_vinterp.py
@@ -0,0 +1,183 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_array_equal
+
+import stratify._bounded_vinterp as bounded_vinterp
+
+
+class Test1D(unittest.TestCase):
+    # 1D cases represent the situation where the vertical coordinate is
+    # 1-dimensional and as such, there are no dimension which the vertical
+    # coordinates vary over other than the 'axis of interpolation'.
+    # The common approach in solving these problems is to reshape and
+    # transpose arrays into the form:
+    # [broadcasting_dims, axis_interpolation, z_varying].
+    # This this will take the form:
+    # [broadcasting_dims, axis_interpolation, 1].
+    def setUp(self):
+        self.data = np.ones((4, 6))
+        self.bounds = self.gen_bounds(0, 7, 1)
+
+    def gen_bounds(self, start, stop, step):
+        bounds = np.vstack([np.arange(start, stop-step, step),
+                            np.arange(start+step, stop, step)])
+        bounds = bounds.transpose((1, 0))
+        return bounds.copy()
+
+    def test_target_half_resolution(self):
+        target_bounds = self.gen_bounds(0, 7, 2)
+        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
+                                          self.data, axis=1)
+        target_data = np.ones((4, 3))*2
+        assert_array_equal(res, target_data)
+
+    def test_target_double_resolution(self):
+        target_bounds = self.gen_bounds(0, 6.5, 0.5)
+        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
+                                          self.data, axis=1)
+        target_data = np.ones((4, 12))*.5
+        assert_array_equal(res, target_data)
+
+    def test_no_broadcasting(self):
+        # In the case of no broadcasting, transposed and reshaped array of form
+        # [broadcasting_dims, axis_interpolation, z_varying] should end up
+        # [1, axis_interpolation, 1].
+        data = self.data[0]
+        target_bounds = self.gen_bounds(0, 7, 2)
+        res = bounded_vinterp.interpolate(target_bounds, self.bounds, data,
+                                          axis=0)
+        target_data = np.ones((3))*2
+        assert_array_equal(res, target_data)
+
+    def test_target_extends_above_source(self):
+        # |-|-|-|-|-|-|-|   - Source
+        # |-|-|-|-|-|-|-|-| - Target
+        source_bounds = self.gen_bounds(0, 7, 1)
+        target_bounds = self.gen_bounds(0, 8, 1)
+        res = bounded_vinterp.interpolate(target_bounds, source_bounds,
+                                          self.data, axis=1)
+        target_data = np.ones((4, 7))
+        target_data[:, -1] = 0
+        assert_array_equal(res, target_data)
+
+    def test_target_extends_below_source(self):
+        #   |-|-|-|-|-|-|-|   - Source
+        # |-|-|-|-|-|-|-|-|   - Target
+        source_bounds = self.gen_bounds(0, 7, 1)
+        target_bounds = self.gen_bounds(-1, 7, 1)
+        res = bounded_vinterp.interpolate(target_bounds, source_bounds,
+                                          self.data, axis=1)
+        target_data = np.ones((4, 7))
+        target_data[:, 0] = 0
+        assert_array_equal(res, target_data)
+
+
+class TestND(unittest.TestCase):
+    # ND cases represent the situation where the vertical coordinate varies
+    # over dimensions other than the axis of interpolation.
+    def setUp(self):
+        self.data = np.ones((2, 6, 4, 3))
+        self.bounds = self.gen_bounds(0, 7, 1)
+
+    def gen_bounds(self, start, stop, step):
+        bounds = np.vstack([np.arange(start, stop-step, step),
+                            np.arange(start+step, stop, step)])
+        bounds = bounds.transpose((1, 0))
+        bounds = bounds[..., None, :].repeat(4, -2)
+        bounds = bounds[..., None, :].repeat(3, -2)
+        return bounds.copy()
+
+    def test_target_half_resolution(self):
+        target_bounds = self.gen_bounds(0, 7, 2)
+        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
+                                          self.data, axis=1)
+
+        target_data = np.ones((2, 3, 4, 3))*2
+        assert_array_equal(res, target_data)
+
+    def test_target_half_resolution_alt_axis(self):
+        # Ensure results as expected with an alternative axis of interpolation.
+        data = self.data.transpose((0, 2, 1, 3))
+        bounds = self.bounds.transpose((1, 0, 2, 3))
+        target_bounds = self.gen_bounds(0, 7, 2)
+        target_bounds = target_bounds.transpose((1, 0, 2, 3))
+
+        res = bounded_vinterp.interpolate(target_bounds, bounds, data, axis=2)
+        target_data = np.ones((2, 4, 3, 3))*2
+        assert_array_equal(res, target_data)
+
+    def test_target_double_resolution(self):
+        target_bounds = self.gen_bounds(0, 6.5, 0.5)
+        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
+                                          self.data, axis=1)
+        target_data = np.ones((2, 12, 4, 3))*.5
+        assert_array_equal(res, target_data)
+
+
+class TestExceptions(unittest.TestCase):
+    def test_mismatch_source_target_level_dimensionality(self):
+        source_bounds = np.zeros((3, 4, 2))
+        target_bounds = np.zeros((4, 2))
+        data = np.zeros((3, 4))
+
+        msg = 'Expecting source and target levels dimensionality'
+        with self.assertRaisesRegexp(ValueError, msg):
+            bounded_vinterp.interpolate(target_bounds, source_bounds, data)
+
+    def test_mismatch_source_target_level_shape(self):
+        # The source and target levels should have identical shape, other than
+        # the axis of interpolation.
+        source_bounds = np.zeros((3, 4, 2))
+        target_bounds = np.zeros((2, 5, 2))
+        data = np.zeros((3, 4))
+
+        msg = 'Expecting the shape of the source and target levels'
+        with self.assertRaisesRegexp(ValueError, msg):
+            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
+                                        axis=0)
+
+    def test_mismatch_between_source_levels_source_data(self):
+        # The source levels should reflect the shape of the data.
+        source_bounds = np.zeros((2, 4, 2))
+        target_bounds = np.zeros((2, 4, 2))
+        data = np.zeros((3, 4))
+
+        msg = 'The provided data is not of compatible shape'
+        with self.assertRaisesRegexp(ValueError, msg):
+            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
+                                        axis=0)
+
+    def test_unexpected_bounds_shape(self):
+        # Expecting bounds of size 2 (that is the upper and lower).
+        # The source levels should reflect the shape of the data.
+        source_bounds = np.zeros((3, 4, 4))
+        target_bounds = np.zeros((4, 4, 4))
+        data = np.zeros((3, 4))
+
+        msg = 'Unexpected source and target bounds shape. shape\[-1\] != 2'
+        with self.assertRaisesRegexp(ValueError, msg):
+            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
+                                        axis=0)
+
+    def test_not_conservative(self):
+        # Where the target does not cover the full extent of the source.
+        # |-|-|-|-|-|-|  - Source
+        #   |-|-|-|-|    - Target
+        def gen_bounds(start, stop, step):
+            bounds = np.vstack([np.arange(start, stop-step, step),
+                                np.arange(start+step, stop, step)])
+            bounds = bounds.transpose((1, 0))
+            return bounds.copy()
+        source_bounds = gen_bounds(0, 7, 1)
+        target_bounds = gen_bounds(1, 6, 1)
+        data = np.ones((4, 6))
+
+        msg = 'Weights calculation yields a less than conservative result.'
+        with self.assertRaisesRegexp(ValueError, msg):
+            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
+                                        axis=1)
+
+
+if __name__ == '__main__':
+    unittest.main()

From 009b38e485dece689e1c52ac2b76448e936774d8 Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Tue, 31 Oct 2017 10:12:38 +0000
Subject: [PATCH 02/10] ENH: Python3 support

---
 stratify/_bounded_vinterp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stratify/_bounded_vinterp.py b/stratify/_bounded_vinterp.py
index 74a5e14..cd61684 100644
--- a/stratify/_bounded_vinterp.py
+++ b/stratify/_bounded_vinterp.py
@@ -94,7 +94,7 @@ def interpolate(z_target, z_src, fz_src, axis=-1):
     # [broadcasting_dims, axis_interpolation, z_varying]
 
     # src_data
-    bdims = range(fz_src.ndim - (z_src.ndim-1))
+    bdims = list(range(fz_src.ndim - (z_src.ndim-1)))
     data_vdims = [ind for ind in range(fz_src.ndim) if ind not in
                   (bdims + [axis])]
     data_transpose = bdims + [axis] + data_vdims

From c87fb50a652ccff39566fa2f57e783ec444dd9b9 Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Tue, 31 Oct 2017 13:02:44 +0000
Subject: [PATCH 03/10] TEST: Flake8 F401 error ignore

---
 stratify/__init__.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/stratify/__init__.py b/stratify/__init__.py
index c1af8c4..281173c 100644
--- a/stratify/__init__.py
+++ b/stratify/__init__.py
@@ -1,11 +1,12 @@
 from __future__ import absolute_import
 
-from ._vinterp import (interpolate, interp_schemes, extrap_schemes,
-                       INTERPOLATE_LINEAR, INTERPOLATE_NEAREST,
+from ._vinterp import (interpolate, interp_schemes,  # noqa: F401
+                       extrap_schemes, INTERPOLATE_LINEAR, INTERPOLATE_NEAREST,
                        EXTRAPOLATE_NAN, EXTRAPOLATE_NEAREST,
-                       EXTRAPOLATE_LINEAR, PyFuncExtrapolator, 
+                       EXTRAPOLATE_LINEAR, PyFuncExtrapolator,
                        PyFuncInterpolator)
-from ._bounded_vinterp import interpolate as interpolate_conservative
+from ._bounded_vinterp import (interpolate as  # noqa: F401
+    interpolate_conservative)
 
 
 __version__ = '0.1a3.dev0'

From a5e9e7c159bbebb430a1ca420a0822c048a042df Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Mon, 6 Nov 2017 16:06:36 +0000
Subject: [PATCH 04/10] MAINT: Review changes

---
 stratify/__init__.py                   |  3 +-
 stratify/_bounded_vinterp.py           |  6 +--
 stratify/_conservative.pyx             |  8 ++--
 stratify/tests/test_bounded_vinterp.py | 61 ++++++++++++++++----------
 4 files changed, 45 insertions(+), 33 deletions(-)

diff --git a/stratify/__init__.py b/stratify/__init__.py
index 281173c..3dbcaae 100644
--- a/stratify/__init__.py
+++ b/stratify/__init__.py
@@ -5,8 +5,7 @@
                        EXTRAPOLATE_NAN, EXTRAPOLATE_NEAREST,
                        EXTRAPOLATE_LINEAR, PyFuncExtrapolator,
                        PyFuncInterpolator)
-from ._bounded_vinterp import (interpolate as  # noqa: F401
-    interpolate_conservative)
+from ._bounded_vinterp import interpolate_conservative  # noqa: F401
 
 
 __version__ = '0.1a3.dev0'
diff --git a/stratify/_bounded_vinterp.py b/stratify/_bounded_vinterp.py
index cd61684..83daaf2 100644
--- a/stratify/_bounded_vinterp.py
+++ b/stratify/_bounded_vinterp.py
@@ -5,9 +5,9 @@
 from ._conservative import conservative_interpolation
 
 
-def interpolate(z_target, z_src, fz_src, axis=-1):
+def interpolate_conservative(z_target, z_src, fz_src, axis=-1):
     """
-    Interface for optimised 1d interpolation across multiple dimensions.
+    1d conservative interpolation across multiple dimensions.
 
     This function provides the ability to perform 1d interpolation on datasets
     with more than one dimension. For instance, this function can be used to
@@ -117,7 +117,7 @@ def interpolate(z_target, z_src, fz_src, axis=-1):
                                              z_src_reshaped.shape[-1]])
     shape = int(np.product(z_target_reshaped.shape[1:-1]))
     z_target_reshaped = z_target_reshaped.reshape(
-            [z_target_reshaped.shape[0], shape, z_target_reshaped.shape[-1]])
+        [z_target_reshaped.shape[0], shape, z_target_reshaped.shape[-1]])
 
     result = conservative_interpolation(
         z_src_reshaped, z_target_reshaped, fz_src_reshaped)
diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
index 7d8069b..351839a 100644
--- a/stratify/_conservative.pyx
+++ b/stratify/_conservative.pyx
@@ -13,8 +13,8 @@ cdef calculate_weights(np.ndarray[np.float64_t, ndim=2] src_point,
     #    weight = overlap / delta
     cdef Py_ssize_t src_ind, tgt_ind
     cdef np.float64_t delta, weight
-    cdef np.ndarray[np.float64_t, ndim = 2] weights
-    cdef np.ndarray[np.float64_t, ndim = 1] src_cell, tgt_cell
+    cdef np.ndarray[np.float64_t, ndim=2] weights
+    cdef np.ndarray[np.float64_t, ndim=1] src_cell, tgt_cell
 
     weights = np.zeros([src_point.shape[0], tgt_point.shape[0]],
                        dtype=np.float64)
@@ -32,8 +32,8 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
                    np.ndarray[np.float64_t, ndim=3] tgt_point,
                    np.ndarray[np.float64_t, ndim=3] src_data):
     cdef Py_ssize_t ind
-    cdef np.ndarray[np.float64_t, ndim = 3] results
-    cdef np.ndarray[np.float64_t, ndim = 2] weights
+    cdef np.ndarray[np.float64_t, ndim=3] results
+    cdef np.ndarray[np.float64_t, ndim=2] weights
     results = np.zeros(
         [src_data.shape[0], tgt_point.shape[0], src_data.shape[2]],
         dtype='float64')
diff --git a/stratify/tests/test_bounded_vinterp.py b/stratify/tests/test_bounded_vinterp.py
index b7d61e9..31c36c6 100644
--- a/stratify/tests/test_bounded_vinterp.py
+++ b/stratify/tests/test_bounded_vinterp.py
@@ -27,15 +27,17 @@ def gen_bounds(self, start, stop, step):
 
     def test_target_half_resolution(self):
         target_bounds = self.gen_bounds(0, 7, 2)
-        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds, self.data,
+                                                       axis=1)
         target_data = np.ones((4, 3))*2
         assert_array_equal(res, target_data)
 
     def test_target_double_resolution(self):
         target_bounds = self.gen_bounds(0, 6.5, 0.5)
-        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds, self.data,
+                                                       axis=1)
         target_data = np.ones((4, 12))*.5
         assert_array_equal(res, target_data)
 
@@ -45,8 +47,9 @@ def test_no_broadcasting(self):
         # [1, axis_interpolation, 1].
         data = self.data[0]
         target_bounds = self.gen_bounds(0, 7, 2)
-        res = bounded_vinterp.interpolate(target_bounds, self.bounds, data,
-                                          axis=0)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds, data,
+                                                       axis=0)
         target_data = np.ones((3))*2
         assert_array_equal(res, target_data)
 
@@ -55,8 +58,9 @@ def test_target_extends_above_source(self):
         # |-|-|-|-|-|-|-|-| - Target
         source_bounds = self.gen_bounds(0, 7, 1)
         target_bounds = self.gen_bounds(0, 8, 1)
-        res = bounded_vinterp.interpolate(target_bounds, source_bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       source_bounds,
+                                                       self.data, axis=1)
         target_data = np.ones((4, 7))
         target_data[:, -1] = 0
         assert_array_equal(res, target_data)
@@ -66,8 +70,9 @@ def test_target_extends_below_source(self):
         # |-|-|-|-|-|-|-|-|   - Target
         source_bounds = self.gen_bounds(0, 7, 1)
         target_bounds = self.gen_bounds(-1, 7, 1)
-        res = bounded_vinterp.interpolate(target_bounds, source_bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       source_bounds,
+                                                       self.data, axis=1)
         target_data = np.ones((4, 7))
         target_data[:, 0] = 0
         assert_array_equal(res, target_data)
@@ -90,8 +95,9 @@ def gen_bounds(self, start, stop, step):
 
     def test_target_half_resolution(self):
         target_bounds = self.gen_bounds(0, 7, 2)
-        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds, self.data,
+                                                       axis=1)
 
         target_data = np.ones((2, 3, 4, 3))*2
         assert_array_equal(res, target_data)
@@ -103,14 +109,16 @@ def test_target_half_resolution_alt_axis(self):
         target_bounds = self.gen_bounds(0, 7, 2)
         target_bounds = target_bounds.transpose((1, 0, 2, 3))
 
-        res = bounded_vinterp.interpolate(target_bounds, bounds, data, axis=2)
+        res = bounded_vinterp.interpolate_conservative(target_bounds, bounds,
+                                                       data, axis=2)
         target_data = np.ones((2, 4, 3, 3))*2
         assert_array_equal(res, target_data)
 
     def test_target_double_resolution(self):
         target_bounds = self.gen_bounds(0, 6.5, 0.5)
-        res = bounded_vinterp.interpolate(target_bounds, self.bounds,
-                                          self.data, axis=1)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds,
+                                                       self.data, axis=1)
         target_data = np.ones((2, 12, 4, 3))*.5
         assert_array_equal(res, target_data)
 
@@ -123,7 +131,8 @@ def test_mismatch_source_target_level_dimensionality(self):
 
         msg = 'Expecting source and target levels dimensionality'
         with self.assertRaisesRegexp(ValueError, msg):
-            bounded_vinterp.interpolate(target_bounds, source_bounds, data)
+            bounded_vinterp.interpolate_conservative(target_bounds,
+                                                     source_bounds, data)
 
     def test_mismatch_source_target_level_shape(self):
         # The source and target levels should have identical shape, other than
@@ -134,8 +143,9 @@ def test_mismatch_source_target_level_shape(self):
 
         msg = 'Expecting the shape of the source and target levels'
         with self.assertRaisesRegexp(ValueError, msg):
-            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
-                                        axis=0)
+            bounded_vinterp.interpolate_conservative(target_bounds,
+                                                     source_bounds, data,
+                                                     axis=0)
 
     def test_mismatch_between_source_levels_source_data(self):
         # The source levels should reflect the shape of the data.
@@ -145,8 +155,9 @@ def test_mismatch_between_source_levels_source_data(self):
 
         msg = 'The provided data is not of compatible shape'
         with self.assertRaisesRegexp(ValueError, msg):
-            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
-                                        axis=0)
+            bounded_vinterp.interpolate_conservative(target_bounds,
+                                                     source_bounds, data,
+                                                     axis=0)
 
     def test_unexpected_bounds_shape(self):
         # Expecting bounds of size 2 (that is the upper and lower).
@@ -157,8 +168,9 @@ def test_unexpected_bounds_shape(self):
 
         msg = 'Unexpected source and target bounds shape. shape\[-1\] != 2'
         with self.assertRaisesRegexp(ValueError, msg):
-            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
-                                        axis=0)
+            bounded_vinterp.interpolate_conservative(target_bounds,
+                                                     source_bounds, data,
+                                                     axis=0)
 
     def test_not_conservative(self):
         # Where the target does not cover the full extent of the source.
@@ -175,8 +187,9 @@ def gen_bounds(start, stop, step):
 
         msg = 'Weights calculation yields a less than conservative result.'
         with self.assertRaisesRegexp(ValueError, msg):
-            bounded_vinterp.interpolate(target_bounds, source_bounds, data,
-                                        axis=1)
+            bounded_vinterp.interpolate_conservative(target_bounds,
+                                                     source_bounds, data,
+                                                     axis=1)
 
 
 if __name__ == '__main__':

From b0a2c5f091caf00b82c670e8602bec6969f56676 Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Mon, 6 Nov 2017 16:19:29 +0000
Subject: [PATCH 05/10] TEST: Extended regexp raised message checks

---
 stratify/tests/test_bounded_vinterp.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/stratify/tests/test_bounded_vinterp.py b/stratify/tests/test_bounded_vinterp.py
index 31c36c6..ff5e178 100644
--- a/stratify/tests/test_bounded_vinterp.py
+++ b/stratify/tests/test_bounded_vinterp.py
@@ -141,7 +141,9 @@ def test_mismatch_source_target_level_shape(self):
         target_bounds = np.zeros((2, 5, 2))
         data = np.zeros((3, 4))
 
-        msg = 'Expecting the shape of the source and target levels'
+        msg = ('Expecting the shape of the source and target levels except '
+               'the axis of interpolation to be identical.  '
+               "\('-', 4, 2\) != \(2, 5, 2\)")
         with self.assertRaisesRegexp(ValueError, msg):
             bounded_vinterp.interpolate_conservative(target_bounds,
                                                      source_bounds, data,
@@ -153,7 +155,8 @@ def test_mismatch_between_source_levels_source_data(self):
         target_bounds = np.zeros((2, 4, 2))
         data = np.zeros((3, 4))
 
-        msg = 'The provided data is not of compatible shape'
+        msg = ('The provided data is not of compatible shape with the '
+               "provided source bounds. \('-', 3, 4\) != \(2, 4\)")
         with self.assertRaisesRegexp(ValueError, msg):
             bounded_vinterp.interpolate_conservative(target_bounds,
                                                      source_bounds, data,

From 44edf42d20d97c4adaf9fbc9ddb87f7acef25b00 Mon Sep 17 00:00:00 2001
From: Carwyn Pelley <cpelley.pub@gmail.com>
Date: Mon, 13 Nov 2017 18:32:42 +0000
Subject: [PATCH 06/10] DOC: Added docstrings

---
 stratify/_conservative.pyx | 76 ++++++++++++++++++++++++++++++++++----
 1 file changed, 69 insertions(+), 7 deletions(-)

diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
index 351839a..5b98e9c 100644
--- a/stratify/_conservative.pyx
+++ b/stratify/_conservative.pyx
@@ -4,13 +4,30 @@ cimport numpy as np
 
 cdef calculate_weights(np.ndarray[np.float64_t, ndim=2] src_point,
                        np.ndarray[np.float64_t, ndim=2] tgt_point):
-    # src_min    src_max
-    #    |----------|                 : Source
-    #       tgt_min     tgt_max
-    #          |------------|         : Target
-    #    |----------|                 : Delta (src_max - src_min)
-    #          |----|                 : Overlap (between src & tgt)
-    #    weight = overlap / delta
+    """
+    Calculate weights for a given point.
+
+    The following visually illustrates the calculation::
+
+        src_min    src_max
+           |----------|                 : Source
+              tgt_min     tgt_max
+                 |------------|         : Target
+           |----------|                 : Delta (src_max - src_min)
+                 |----|                 : Overlap (between src & tgt)
+           weight = overlap / delta
+
+    Parameters
+    ----------
+    src_point (2d double array) - Source point (at a specific location).
+    tgt_point (2d double array) - Target point (at a specific location).
+
+    Returns
+    -------
+    2d double array - Weights corresponding to shape [src_point.shape[0],
+        tgt_point.shape[0]].
+
+    """
     cdef Py_ssize_t src_ind, tgt_ind
     cdef np.float64_t delta, weight
     cdef np.ndarray[np.float64_t, ndim=2] weights
@@ -31,6 +48,28 @@ cdef calculate_weights(np.ndarray[np.float64_t, ndim=2] src_point,
 cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
                    np.ndarray[np.float64_t, ndim=3] tgt_point,
                    np.ndarray[np.float64_t, ndim=3] src_data):
+    """
+    Perform conservative interpolation.
+
+    Conservation interpolation on of a dataset between a provided source
+    coordinate and a target coordinate.
+
+    Parameters
+    ----------
+    src_points (3d double array) - Source coordinate, taking the form
+        [axis_interpolation, z_varying, 2].
+    tgt_points (3d double array) - Target coordinate, taking the form
+        [axis_interpolation, z_varying, 2].
+    src_data (3d double array) - The source data, the phenonenon data to be
+        interpolated from ``src_points`` to ``tgt_points``.  Taking the form
+        [broadcasting_dims, axis_interpolation, z_varying].
+
+    Returns
+    -------
+    3d double array - Interpolated result over target levels (``tgt_points``).
+        Taking the form [broadcasting_dims, axis_interpolation, z_varying].
+
+    """
     cdef Py_ssize_t ind
     cdef np.ndarray[np.float64_t, ndim=3] results
     cdef np.ndarray[np.float64_t, ndim=2] weights
@@ -50,6 +89,29 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
 
 
 def conservative_interpolation(src_points, tgt_points, src_data):
+    """
+    Perform conservative interpolation.
+
+    Conservation interpolation of a dataset between a provided source
+    coordinate and a target coordinate.  All inputs are recast to 64-bit float
+    arrays before calculation.
+
+    Parameters
+    ----------
+    src_points (3d array) - Source coordinate, taking the form
+        [axis_interpolation, z_varying, 2].
+    tgt_points (3d array) - Target coordinate, taking the form
+        [axis_interpolation, z_varying, 2].
+    src_data (3d array) - The source data, the phenonenon data to be
+        interpolated from ``src_points`` to ``tgt_points``.  Taking the form
+        [broadcasting_dims, axis_interpolation, z_varying].
+
+    Returns
+    -------
+    3d double array - Interpolated result over target levels (``tgt_points``).
+        Taking the form [broadcasting_dims, axis_interpolation, z_varying].
+
+    """
     return apply_weights(src_points.astype('float64'),
                          tgt_points.astype('float64'),
                          src_data.astype('float64'))

From 20eb3f5fa4822ca5eb4e40d690aa7cd92e4f87f0 Mon Sep 17 00:00:00 2001
From: cpelley <cpelley.pub@gmail.com>
Date: Tue, 16 Jan 2018 18:14:31 +0000
Subject: [PATCH 07/10] TEST: Added non-equally spaced coordinate test

---
 stratify/tests/test_bounded_vinterp.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/stratify/tests/test_bounded_vinterp.py b/stratify/tests/test_bounded_vinterp.py
index ff5e178..e143c1e 100644
--- a/stratify/tests/test_bounded_vinterp.py
+++ b/stratify/tests/test_bounded_vinterp.py
@@ -65,6 +65,20 @@ def test_target_extends_above_source(self):
         target_data[:, -1] = 0
         assert_array_equal(res, target_data)
 
+    def test_target_extends_above_source_non_equally_spaced_coords(self):
+        # |--|--|-------||  - Source
+        # |-|-|-|-|-|-|-|-| - Target
+        source_bounds = np.array([[0, 1.5], [1.5, 2], [2, 6], [6, 6.5]])
+        target_bounds = self.gen_bounds(0, 8, 1)
+        data = np.ones((4, 4))
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       source_bounds, data,
+                                                       axis=1)
+        target_data = np.array(
+            [1/1.5, 1 + ((1/3.)/1), 0.25, 0.25, 0.25, 0.25, 1.])[None]
+        target_data = np.repeat(target_data, 4, 0)
+        assert_array_equal(res, target_data)
+
     def test_target_extends_below_source(self):
         #   |-|-|-|-|-|-|-|   - Source
         # |-|-|-|-|-|-|-|-|   - Target

From 54bb77a9efc9fe3748cdfca0fd31f854d3f82c22 Mon Sep 17 00:00:00 2001
From: cpelley <cpelley.pub@gmail.com>
Date: Tue, 16 Jan 2018 18:27:34 +0000
Subject: [PATCH 08/10] ENH: Return nan values where no source cells contribute
 to a target cell.

TODO: Update tests
---
 stratify/_conservative.pyx | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
index 5b98e9c..46c1038 100644
--- a/stratify/_conservative.pyx
+++ b/stratify/_conservative.pyx
@@ -51,8 +51,9 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
     """
     Perform conservative interpolation.
 
-    Conservation interpolation on of a dataset between a provided source
-    coordinate and a target coordinate.
+    Conservation interpolation of a dataset between a provided source
+    coordinate and a target coordinate.  Where no source cells contribute to a
+    target cell, a np.nan value is returned.
 
     Parameters
     ----------
@@ -60,7 +61,7 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
         [axis_interpolation, z_varying, 2].
     tgt_points (3d double array) - Target coordinate, taking the form
         [axis_interpolation, z_varying, 2].
-    src_data (3d double array) - The source data, the phenonenon data to be
+    src_data (3d double array) - The source data, the phenomenon data to be
         interpolated from ``src_points`` to ``tgt_points``.  Taking the form
         [broadcasting_dims, axis_interpolation, z_varying].
 
@@ -85,6 +86,8 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
             raise ValueError(msg)
         results[..., ind] = (
             weights * src_data[..., ind][..., None]).sum(axis=1)
+    # Return np.nan for those target cells where no source contributes.
+    results[:, weights.sum(axis=0) == 0, :] = np.nan
     return results
 
 

From 58098b81c3c2b549ce37f94e181b6fc6db25ab23 Mon Sep 17 00:00:00 2001
From: cpelley <cpelley.pub@gmail.com>
Date: Wed, 17 Jan 2018 12:22:22 +0000
Subject: [PATCH 09/10] ENH: Return nan values in tgt cells where a src nan
 cells contribute

---
 stratify/_bounded_vinterp.py           |  4 +++-
 stratify/_conservative.pyx             | 15 +++++++++++----
 stratify/tests/test_bounded_vinterp.py | 19 +++++++++++++++++--
 3 files changed, 31 insertions(+), 7 deletions(-)

diff --git a/stratify/_bounded_vinterp.py b/stratify/_bounded_vinterp.py
index 83daaf2..b2c865c 100644
--- a/stratify/_bounded_vinterp.py
+++ b/stratify/_bounded_vinterp.py
@@ -49,8 +49,10 @@ def interpolate_conservative(z_target, z_src, fz_src, axis=-1):
 
     Note
     ----
-    Support for 1D z_target and corresponding ND z_src will be provided in
+    - Support for 1D z_target and corresponding ND z_src will be provided in
     future as driven by user requirement.
+    - Those cells, where 'nan' values in the source data contribute, a 'nan'
+      value is returned.
 
     """
     if z_src.ndim != z_target.ndim:
diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
index 46c1038..ebe6786 100644
--- a/stratify/_conservative.pyx
+++ b/stratify/_conservative.pyx
@@ -72,8 +72,9 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
 
     """
     cdef Py_ssize_t ind
-    cdef np.ndarray[np.float64_t, ndim=3] results
+    cdef np.ndarray[np.float64_t, ndim=3] results, weighted_contrib
     cdef np.ndarray[np.float64_t, ndim=2] weights
+    cdef np.ndarray[np.int64_t, ndim=2] nan_src_contrib
     results = np.zeros(
         [src_data.shape[0], tgt_point.shape[0], src_data.shape[2]],
         dtype='float64')
@@ -84,10 +85,16 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
             msg = ('Weights calculation yields a less than conservative '
                    'result.  Aborting.')
             raise ValueError(msg)
+        weighted_contrib = weights * src_data[..., ind][..., None]
         results[..., ind] = (
-            weights * src_data[..., ind][..., None]).sum(axis=1)
-    # Return np.nan for those target cells where no source contributes.
-    results[:, weights.sum(axis=0) == 0, :] = np.nan
+            np.nansum(weighted_contrib, axis=1))
+        # Return nan values for those target cells, where there is any
+        # contribution of nan data from the source data.
+        nan_src_contrib = ((weights > 0) * np.isnan(weighted_contrib)).sum(axis=1)
+        results[..., ind][nan_src_contrib>0] = np.nan
+
+        # Return np.nan for those target cells where no source contributes.
+        results[:, weights.sum(axis=0) == 0, ind] = np.nan
     return results
 
 
diff --git a/stratify/tests/test_bounded_vinterp.py b/stratify/tests/test_bounded_vinterp.py
index e143c1e..58c48c7 100644
--- a/stratify/tests/test_bounded_vinterp.py
+++ b/stratify/tests/test_bounded_vinterp.py
@@ -53,6 +53,21 @@ def test_no_broadcasting(self):
         target_data = np.ones((3))*2
         assert_array_equal(res, target_data)
 
+    def test_source_with_nans(self):
+        # In the case of no broadcasting, transposed and reshaped array of form
+        # [broadcasting_dims, axis_interpolation, z_varying] should end up
+        # [1, axis_interpolation, 1].
+        data = self.data[0]
+        data[0] = data[-2:] = np.nan
+        target_bounds = self.gen_bounds(0, 6.5, 0.5)
+        res = bounded_vinterp.interpolate_conservative(target_bounds,
+                                                       self.bounds, data,
+                                                       axis=0)
+        target_data = np.ones((12))*.5
+        target_data[:2] = np.nan
+        target_data[-4:] = np.nan
+        assert_array_equal(res, target_data)
+
     def test_target_extends_above_source(self):
         # |-|-|-|-|-|-|-|   - Source
         # |-|-|-|-|-|-|-|-| - Target
@@ -62,7 +77,7 @@ def test_target_extends_above_source(self):
                                                        source_bounds,
                                                        self.data, axis=1)
         target_data = np.ones((4, 7))
-        target_data[:, -1] = 0
+        target_data[:, -1] = np.nan
         assert_array_equal(res, target_data)
 
     def test_target_extends_above_source_non_equally_spaced_coords(self):
@@ -88,7 +103,7 @@ def test_target_extends_below_source(self):
                                                        source_bounds,
                                                        self.data, axis=1)
         target_data = np.ones((4, 7))
-        target_data[:, 0] = 0
+        target_data[:, 0] = np.nan
         assert_array_equal(res, target_data)
 
 

From bab0fdd533092fe94a20da3f6c90881ef9274f1f Mon Sep 17 00:00:00 2001
From: cpelley <cpelley.pub@gmail.com>
Date: Wed, 17 Jan 2018 12:44:40 +0000
Subject: [PATCH 10/10] MAINT: Slight simplification

---
 stratify/_conservative.pyx | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/stratify/_conservative.pyx b/stratify/_conservative.pyx
index ebe6786..ab8e0e0 100644
--- a/stratify/_conservative.pyx
+++ b/stratify/_conservative.pyx
@@ -74,7 +74,6 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
     cdef Py_ssize_t ind
     cdef np.ndarray[np.float64_t, ndim=3] results, weighted_contrib
     cdef np.ndarray[np.float64_t, ndim=2] weights
-    cdef np.ndarray[np.int64_t, ndim=2] nan_src_contrib
     results = np.zeros(
         [src_data.shape[0], tgt_point.shape[0], src_data.shape[2]],
         dtype='float64')
@@ -90,8 +89,8 @@ cdef apply_weights(np.ndarray[np.float64_t, ndim=3] src_point,
             np.nansum(weighted_contrib, axis=1))
         # Return nan values for those target cells, where there is any
         # contribution of nan data from the source data.
-        nan_src_contrib = ((weights > 0) * np.isnan(weighted_contrib)).sum(axis=1)
-        results[..., ind][nan_src_contrib>0] = np.nan
+        results[..., ind][
+            ((weights > 0) * np.isnan(weighted_contrib)).any(axis=1)] = np.nan
 
         # Return np.nan for those target cells where no source contributes.
         results[:, weights.sum(axis=0) == 0, ind] = np.nan