Index: /trunk/anuga_core/source/anuga/fit_interpolate/interpolate2d.py
===================================================================
--- /trunk/anuga_core/source/anuga/fit_interpolate/interpolate2d.py	(revision 8967)
+++ /trunk/anuga_core/source/anuga/fit_interpolate/interpolate2d.py	(revision 8967)
@@ -0,0 +1,302 @@
+"""Bilinear and piecewise constant interpolation.
+
+* Interpolate from a regular numpy array to arbitrary points.
+* NaN's are taken care of sensibly.
+* Algorithm works both for cartesian (x, y) and geographic (lon, lat)
+  coordinates.
+* Guarantees interpolated values within extrema of grid points, i.e.
+  no overshoot.
+
+See docstring for details on usage and the mathematical derivation.
+
+To test, run the accompanying script test_interpolate2d.py
+
+Author: Ole Nielsen 2011
+
+"""
+__author__ = 'Ole Nielsen <ole.moller.nielsen@gmail.com>'
+__revision__ = '$Format:%H$'
+__date__ = '01/11/2011'
+__license__ = 'GPL v3'
+__copyright__ = 'Copyright 2012, Australia Indonesia Facility for '
+__copyright__ += 'Disaster Reduction'
+
+import numpy
+from anuga.anuga_exceptions import ANUGAError
+
+
+def interpolate2d(x, y, Z, points, mode='linear', bounds_error=False):
+    """Fundamental 2D interpolation routine
+
+    Input
+        x: 1D array of x-coordinates of the mesh on which to interpolate
+        y: 1D array of y-coordinates of the mesh on which to interpolate
+        Z: 2D array of values for each x, y pair
+        points: Nx2 array of coordinates where interpolated values are sought
+        mode: Determines the interpolation order. Options are
+            * 'constant' - piecewise constant nearest neighbour interpolation
+            * 'linear' - bilinear interpolation using the four
+                         nearest neighbours (default)
+        bounds_error: Boolean flag. If True (default) a BoundsError exception
+            will be raised when interpolated values are requested
+            outside the domain of the input data. If False, nan
+            is returned for those values
+
+    Output
+        1D array of length N (same length as points) with interpolated values
+
+    Notes
+        Input coordinates x and y are assumed to be monotonically increasing,
+        but need not be equidistantly spaced.
+
+        Z is assumed to have dimension M x N, where M = len(x) and N = len(y).
+        In other words it is assumed that the x values follow the first
+        (vertical) axis downwards and y values the second (horizontal) axis
+        from left to right.
+
+        If this routine is to be used for interpolation of raster grids where
+        data is typically organised with longitudes (x) going from left to
+        right and latitudes (y) from left to right then user
+        interpolate_raster in this module
+
+    Derivation
+
+        Bilinear interpolation is based on the standard 1D linear interpolation
+        formula:
+
+        Given points (x0, y0) and (x1, y1) and a value of x where
+        x0 <= x <= x1, the linearly interpolated value y at x is given as
+
+        alpha*(y1-y0) + y0
+
+        or
+
+        alpha*y1 + (1-alpha)*y0                (1)
+
+        where alpha = (x-x0)/(x1-x0)           (1a)
+
+
+        2D bilinear interpolation aims at obtaining an interpolated value z at
+        a point (x,y) which lies inside a square formed by points (x0, y0),
+        (x1, y0), (x0, y1) and (x1, y1) for which values z00, z10, z01 and
+        z11 are known.
+
+        This is obtained by applying equation (1) twice  in the
+        x-direction to obtain interpolated points q0 and q1 for (x, y0) and
+        (x, y1), respectively.
+
+        q0 = alpha*z10 + (1-alpha)*z00         (2)
+
+        and
+
+        q1 = alpha*z11 + (1-alpha)*z01         (3)
+
+
+        Then using equation (1) in the y-direction on the results from (2)
+        and (3)
+
+        z = beta*q1 + (1-beta)*q0              (4)
+
+        where beta = (y-y0)/(y1-y0)            (4a)
+
+
+        Substituting (2) and (3) into (4) yields
+
+        z = alpha*beta*z11 + beta*z01 - alpha*beta*z01 +
+            alpha*z10 + z00 - alpha*z00 - alpha*beta*z10 - beta*z00 +
+            alpha*beta*z00
+          = alpha*beta*(z11 - z01 - z10 + z00) +
+            alpha*(z10 - z00) + beta*(z01 - z00) + z00
+
+        which can be further simplified to
+
+        z = alpha*beta*(z11 - dx - dy - z00) + alpha*dx + beta*dy + z00  (5)
+
+        where
+        dx = z10 - z00
+        dy = z01 - z00
+
+        Equation (5) is what is implemented in the function interpolate2d.
+
+
+        Piecewise constant interpolation can be implemented using the same
+        coefficients (1a) and (4a) that are used for bilinear interpolation
+        as they are a measure of the relative distance to the left and lower
+        neigbours. A value of 0 will pick the left or lower bound whereas a
+        value of 1 will pick the right or higher bound. Hence z can be
+        assigned to its nearest neigbour as follows
+
+            | z00   alpha < 0.5 and beta < 0.5    # lower left corner
+            |
+            | z10   alpha >= 0.5 and beta < 0.5   # lower right corner
+        z = |
+            | z01   alpha < 0.5 and beta >= 0.5   # upper left corner
+            |
+            | z11   alpha >= 0.5 and beta >= 0.5  # upper right corner
+
+    """
+
+    # Check inputs and provide xi, eta as x and y coordinates from
+    # points vector
+    x, y, Z, xi, eta = check_inputs(x, y, Z, points, mode, bounds_error)
+
+    # Identify elements that are outside interpolation domain or NaN
+    outside = (xi < x[0]) + (eta < y[0]) + (xi > x[-1]) + (eta > y[-1])
+    outside += numpy.isnan(xi) + numpy.isnan(eta)
+
+    # Restrict interpolation points to those that are inside the grid
+    inside = -outside  # Invert boolean array to find elements inside
+    xi = xi[inside]
+    eta = eta[inside]
+
+    # Find upper neighbours for each interpolation point
+    # ('left' means first occurrence)
+    idx = numpy.searchsorted(x, xi, side='left')
+    idy = numpy.searchsorted(y, eta, side='left')
+
+    # Get the four neighbours for each interpolation point
+    x0 = x[idx - 1]  # Left
+    x1 = x[idx]      # Right
+    y0 = y[idy - 1]  # Lower
+    y1 = y[idy]      # Upper
+
+    # And the corresponding four grid values
+    z00 = Z[idx - 1, idy - 1]
+    z01 = Z[idx - 1, idy]
+    z10 = Z[idx, idy - 1]
+    z11 = Z[idx, idy]
+
+    # Coefficients for weighting between lower and upper bounds
+    oldset = numpy.seterr(invalid='ignore')  # Suppress zero division warning
+    alpha = (xi - x0) / (x1 - x0)
+    beta = (eta - y0) / (y1 - y0)
+    numpy.seterr(**oldset)  # Restore warnings
+
+    if mode == 'linear':
+        # Bilinear interpolation formula as per equation (5) above
+        dx = z10 - z00
+        dy = z01 - z00
+        z = z00 + alpha * dx + beta * dy + alpha * beta * (z11 - dx - dy - z00)
+    else:
+        # Piecewise constant (as verified in input_check)
+
+        # Set up masks for the quadrants
+        left = alpha < 0.5
+        right = -left
+        lower = beta < 0.5
+        upper = -lower
+
+        lower_left = lower * left
+        lower_right = lower * right
+        upper_left = upper * left
+
+        # Initialise result array with all elements set to upper right
+        z = z11
+
+        # Then set the other quadrants
+        z[lower_left] = z00[lower_left]
+        z[lower_right] = z10[lower_right]
+        z[upper_left] = z01[upper_left]
+
+    # Populate result with interpolated values for points inside domain
+    # and NaN for values outside
+    r = numpy.zeros(len(points))
+    r[inside] = z
+    r[outside] = numpy.nan
+
+    # Return interpolated points
+    return r
+
+
+#------------------------
+# Auxiliary functionality
+#------------------------
+def check_inputs(x, y, Z, points, mode, bounds_error):
+    """Check inputs for interpolate2d function
+    """
+
+    msg = ('Only mode "linear" and "constant" are implemented. '
+           'I got "%s"' % mode)
+    if mode not in ['linear', 'constant']:
+        raise ANUGAError(msg)
+
+    x = numpy.array(x)
+
+    try:
+        y = numpy.array(y)
+    except Exception, e:
+        msg = ('Input vector y could not be converted to numpy array: '
+               '%s' % str(e))
+        raise Exception(msg)
+
+    msg = ('Input vector x must be monotonically increasing. I got '
+           'min(x) == %.15f, but x[0] == %.15f' % (min(x), x[0]))
+    if not min(x) == x[0]:
+        raise ANUGAError(msg)
+
+    msg = ('Input vector y must be monotoneously increasing. '
+           'I got min(y) == %.15f, but y[0] == %.15f' % (min(y), y[0]))
+    if not min(y) == y[0]:
+        raise ANUGAError(msg)
+
+    msg = ('Input vector x must be monotoneously increasing. I got '
+           'max(x) == %.15f, but x[-1] == %.15f' % (max(x), x[-1]))
+    if not max(x) == x[-1]:
+        raise ANUGAError(msg)
+
+    msg = ('Input vector y must be monotoneously increasing. I got '
+           'max(y) == %.15f, but y[-1] == %.15f' % (max(y), y[-1]))
+    if not max(y) == y[-1]:
+        raise ANUGAError(msg)
+
+    try:
+        Z = numpy.array(Z)
+        m, n = Z.shape
+    except Exception, e:
+        msg = 'Z must be a 2D numpy array: %s' % str(e)
+        raise Exception(msg)
+
+    Nx = len(x)
+    Ny = len(y)
+    msg = ('Input array Z must have dimensions %i x %i corresponding to the '
+           'lengths of the input coordinates x and y. However, '
+           'Z has dimensions %i x %i.' % (Nx, Ny, m, n))
+    if not (Nx == m and Ny == n):
+        raise ANUGAError(msg)
+
+    # Get interpolation points
+    points = numpy.array(points)
+    xi = points[:, 0]
+    eta = points[:, 1]
+
+    if bounds_error:
+        xi0 = min(xi)
+        xi1 = max(xi)
+        eta0 = min(eta)
+        eta1 = max(eta)
+
+        msg = ('Interpolation point xi=%f was less than the smallest '
+               'value in domain (x=%f) and bounds_error was requested.'
+               % (xi0, x[0]))
+        if xi0 < x[0]:
+            raise BoundsError(msg)
+
+        msg = ('Interpolation point xi=%f was greater than the largest '
+               'value in domain (x=%f) and bounds_error was requested.'
+               % (xi1, x[-1]))
+        if xi1 > x[-1]:
+            raise BoundsError(msg)
+
+        msg = ('Interpolation point eta=%f was less than the smallest '
+               'value in domain (y=%f) and bounds_error was requested.'
+               % (eta0, y[0]))
+        if eta0 < y[0]:
+            raise BoundsError(msg)
+
+        msg = ('Interpolation point eta=%f was greater than the largest '
+               'value in domain (y=%f) and bounds_error was requested.'
+               % (eta1, y[-1]))
+        if eta1 > y[-1]:
+            raise BoundsError(msg)
+
+    return x, y, Z, xi, eta
Index: /trunk/anuga_core/source/anuga/fit_interpolate/test_interpolate2d.py
===================================================================
--- /trunk/anuga_core/source/anuga/fit_interpolate/test_interpolate2d.py	(revision 8967)
+++ /trunk/anuga_core/source/anuga/fit_interpolate/test_interpolate2d.py	(revision 8967)
@@ -0,0 +1,466 @@
+"""Test suite for interpolate2d.py
+
+Ole Nielsen, 2011
+"""
+
+__author__ = 'Ole Nielsen <ole.moller.nielsen@gmail.com>'
+__revision__ = '$Format:%H$'
+__date__ = '01/11/2011'
+__license__ = 'GPL v3'
+__copyright__ = 'Copyright 2012, Australia Indonesia Facility for '
+__copyright__ += 'Disaster Reduction'
+
+import unittest
+import numpy
+
+from interpolate2d import interpolate2d
+
+
+# Auxiliary functions for testing
+def nanallclose(x, y, rtol=1.0e-5, atol=1.0e-8):
+    """Numpy allclose function which allows NaN
+
+    Args:
+        * x, y: Either scalars or numpy arrays
+
+    Returns:
+        * True or False
+
+    Note:
+        Returns True if all non-nan elements pass.
+    """
+
+    xn = numpy.isnan(x)
+    yn = numpy.isnan(y)
+    if numpy.any(xn != yn):
+        # Presence of NaNs is not the same in x and y
+        return False
+
+    if numpy.all(xn):
+        # Everything is NaN.
+        # This will also take care of x and y being NaN scalars
+        return True
+
+    # Filter NaN's out
+    if numpy.any(xn):
+        x = x[-xn]
+        y = y[-yn]
+
+    # Compare non NaN's and return
+    return numpy.allclose(x, y, rtol=rtol, atol=atol)
+
+
+def axes2points(x, y):
+    """Generate all combinations of grid point coordinates from x and y axes
+
+    Args:
+        * x: x coordinates (array)
+        * y: y coordinates (array)
+
+    Returns:
+        * P: Nx2 array consisting of coordinates for all
+             grid points defined by x and y axes. The x coordinate
+             will vary the fastest to match the way 2D numpy
+             arrays are laid out by default ('C' order). That way,
+             the x and y coordinates will match a corresponding
+             2D array A when flattened (A.flat[:] or A.reshape(-1))
+
+    Note:
+        Example
+
+        x = [1, 2, 3]
+        y = [10, 20]
+
+        P = [[1, 10],
+             [2, 10],
+             [3, 10],
+             [1, 20],
+             [2, 20],
+             [3, 20]]
+    """
+
+    # Reverse y coordinates to have them start at bottom of array
+    y = numpy.flipud(y)
+
+    # Repeat x coordinates for each y (fastest varying)
+    X = numpy.kron(numpy.ones(len(y)), x)
+
+    # Repeat y coordinates for each x (slowest varying)
+    Y = numpy.kron(y, numpy.ones(len(x)))
+
+    # Check
+    N = len(X)
+    assert len(Y) == N
+
+    # Create Nx2 array of x and y coordinates
+    X = numpy.reshape(X, (N, 1))
+    Y = numpy.reshape(Y, (N, 1))
+    P = numpy.concatenate((X, Y), axis=1)
+
+    # Return
+    return P
+
+
+def linear_function(x, y):
+    """Auxiliary function for use with interpolation test
+    """
+
+    return x + y / 2.0
+
+
+class Test_interpolate(unittest.TestCase):
+
+    def test_linear_interpolation_basic(self):
+        """Interpolation library works for linear function - basic test
+        """
+
+        # Define pixel centers along each direction
+        x = [1.0, 2.0, 4.0]
+        y = [5.0, 9.0]
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define values for each x, y pair as a linear function
+        for i in range(len(x)):
+            for j in range(len(y)):
+                A[i, j] = linear_function(x[i], y[j])
+
+        # Test first that original points are reproduced correctly
+        for i, xi in enumerate(x):
+            for j, eta in enumerate(y):
+                val = interpolate2d(x, y, A, [(xi, eta)], mode='linear')[0]
+                ref = linear_function(xi, eta)
+                assert numpy.allclose(val, ref, rtol=1e-12, atol=1e-12)
+
+        # Then test that genuinly interpolated points are correct
+        xis = numpy.linspace(x[0], x[-1], 10)
+        etas = numpy.linspace(y[0], y[-1], 10)
+        points = axes2points(xis, etas)
+
+        vals = interpolate2d(x, y, A, points, mode='linear')
+        refs = linear_function(points[:, 0], points[:, 1])
+        assert numpy.allclose(vals, refs, rtol=1e-12, atol=1e-12)
+
+    def test_constant_interpolation_basic(self):
+        """Interpolation library works for piecewise constant function
+        """
+
+        # Define pixel centers along each direction
+        x = numpy.array([1.0, 2.0, 4.0])
+        y = numpy.array([5.0, 9.0])
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define values for each x, y pair as a linear function
+        for i in range(len(x)):
+            for j in range(len(y)):
+                A[i, j] = linear_function(x[i], y[j])
+
+        # Then test that interpolated points are always assigned value of
+        # closest neighbour
+        xis = numpy.linspace(x[0], x[-1], 10)
+        etas = numpy.linspace(y[0], y[-1], 10)
+        points = axes2points(xis, etas)
+
+        vals = interpolate2d(x, y, A, points, mode='constant')
+
+        # Find upper neighbours for each interpolation point
+        xi = points[:, 0]
+        eta = points[:, 1]
+        idx = numpy.searchsorted(x, xi, side='left')
+        idy = numpy.searchsorted(y, eta, side='left')
+
+        # Get the four neighbours for each interpolation point
+        x0 = x[idx - 1]
+        x1 = x[idx]
+        y0 = y[idy - 1]
+        y1 = y[idy]
+
+        z00 = A[idx - 1, idy - 1]
+        z01 = A[idx - 1, idy]
+        z10 = A[idx, idy - 1]
+        z11 = A[idx, idy]
+
+        # Location coefficients
+        alpha = (xi - x0) / (x1 - x0)
+        beta = (eta - y0) / (y1 - y0)
+
+        refs = numpy.zeros(len(vals))
+        for i in range(len(refs)):
+            if alpha[i] < 0.5 and beta[i] < 0.5:
+                refs[i] = z00[i]
+
+            if alpha[i] >= 0.5 and beta[i] < 0.5:
+                refs[i] = z10[i]
+
+            if alpha[i] < 0.5 and beta[i] >= 0.5:
+                refs[i] = z01[i]
+
+            if alpha[i] >= 0.5 and beta[i] >= 0.5:
+                refs[i] = z11[i]
+
+        assert numpy.allclose(vals, refs, rtol=1e-12, atol=1e-12)
+
+    def test_linear_interpolation_range(self):
+        """Interpolation library works for linear function - a range of cases
+        """
+
+        for x in [[1.0, 2.0, 4.0], [-20, -19, 0], numpy.arange(200) + 1000]:
+            for y in [[5.0, 9.0], [100, 200, 10000]]:
+
+                # Define ny by nx array with corresponding values
+                A = numpy.zeros((len(x), len(y)))
+
+                # Define values for each x, y pair as a linear function
+                for i in range(len(x)):
+                    for j in range(len(y)):
+                        A[i, j] = linear_function(x[i], y[j])
+
+                # Test that linearly interpolated points are correct
+                xis = numpy.linspace(x[0], x[-1], 100)
+                etas = numpy.linspace(y[0], y[-1], 100)
+                points = axes2points(xis, etas)
+
+                vals = interpolate2d(x, y, A, points, mode='linear')
+                refs = linear_function(points[:, 0], points[:, 1])
+                assert numpy.allclose(vals, refs, rtol=1e-12, atol=1e-12)
+
+    def test_linear_interpolation_nan_points(self):
+        """Interpolation library works with interpolation points being NaN
+
+        This is was the reason for bug reported in:
+        https://github.com/AIFDR/riab/issues/155
+        """
+
+        # Define pixel centers along each direction
+        x = [1.0, 2.0, 4.0]
+        y = [5.0, 9.0]
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define values for each x, y pair as a linear function
+        for i in range(len(x)):
+            for j in range(len(y)):
+                A[i, j] = linear_function(x[i], y[j])
+
+        # Then test that interpolated points can contain NaN
+        xis = numpy.linspace(x[0], x[-1], 10)
+        etas = numpy.linspace(y[0], y[-1], 10)
+        xis[6:7] = numpy.nan
+        etas[3] = numpy.nan
+        points = axes2points(xis, etas)
+
+        vals = interpolate2d(x, y, A, points, mode='linear')
+        refs = linear_function(points[:, 0], points[:, 1])
+        assert nanallclose(vals, refs, rtol=1e-12, atol=1e-12)
+
+    def test_linear_interpolation_nan_array(self):
+        """Interpolation library works (linear mode) with grid points being NaN
+        """
+
+        # Define pixel centers along each direction
+        x = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
+        y = [4.0, 5.0, 7.0, 9.0, 11.0, 13.0]
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define values for each x, y pair as a linear function
+        for i in range(len(x)):
+            for j in range(len(y)):
+                A[i, j] = linear_function(x[i], y[j])
+        A[2, 3] = numpy.nan  # (x=2.0, y=9.0): NaN
+
+        # Then test that interpolated points can contain NaN
+        xis = numpy.linspace(x[0], x[-1], 12)
+        etas = numpy.linspace(y[0], y[-1], 10)
+        points = axes2points(xis, etas)
+
+        vals = interpolate2d(x, y, A, points, mode='linear')
+        refs = linear_function(points[:, 0], points[:, 1])
+
+        # Set reference result with expected NaNs and compare
+        for i, (xi, eta) in enumerate(points):
+            if (1.0 < xi <= 3.0) and (7.0 < eta <= 11.0):
+                refs[i] = numpy.nan
+
+        assert nanallclose(vals, refs, rtol=1e-12, atol=1e-12)
+
+    def test_interpolation_random_array_and_nan(self):
+        """Interpolation library (constant and linear) works with NaN
+        """
+
+        # Define pixel centers along each direction
+        x = numpy.arange(20) * 1.0
+        y = numpy.arange(25) * 1.0
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define arbitrary values for each x, y pair
+        numpy.random.seed(17)
+        A = numpy.random.random((len(x), len(y))) * 10
+
+        # Create islands of NaN
+        A[5, 13] = numpy.nan
+        A[6, 14] = A[6, 18] = numpy.nan
+        A[7, 14:18] = numpy.nan
+        A[8, 13:18] = numpy.nan
+        A[9, 12:19] = numpy.nan
+        A[10, 14:17] = numpy.nan
+        A[11, 15] = numpy.nan
+
+        A[15, 5:6] = numpy.nan
+
+        # Creat interpolation points
+        xis = numpy.linspace(x[0], x[-1], 39)   # Hit all mid points
+        etas = numpy.linspace(y[0], y[-1], 73)  # Hit thirds
+        points = axes2points(xis, etas)
+
+        for mode in ['linear', 'constant']:
+            vals = interpolate2d(x, y, A, points, mode=mode)
+
+            # Calculate reference result with expected NaNs and compare
+            i = j = 0
+            for k, (xi, eta) in enumerate(points):
+
+                # Find indices of nearest higher value in x and y
+                i = numpy.searchsorted(x, xi)
+                j = numpy.searchsorted(y, eta)
+
+                if i > 0 and j > 0:
+
+                    # Get four neigbours
+                    A00 = A[i - 1, j - 1]
+                    A01 = A[i - 1, j]
+                    A10 = A[i, j - 1]
+                    A11 = A[i, j]
+
+                    if numpy.allclose(xi, x[i]):
+                        alpha = 1.0
+                    else:
+                        alpha = 0.5
+
+                    if numpy.allclose(eta, y[j]):
+                        beta = 1.0
+                    else:
+                        beta = eta - y[j - 1]
+
+                    if mode == 'linear':
+                        if numpy.any(numpy.isnan([A00, A01, A10, A11])):
+                            ref = numpy.nan
+                        else:
+                            ref = (A00 * (1 - alpha) * (1 - beta) +
+                                   A01 * (1 - alpha) * beta +
+                                   A10 * alpha * (1 - beta) +
+                                   A11 * alpha * beta)
+                    elif mode == 'constant':
+                        assert alpha >= 0.5  # Only case in this test
+
+                        if beta < 0.5:
+                            ref = A10
+                        else:
+                            ref = A11
+                    else:
+                        msg = 'Unknown mode: %s' % mode
+                        raise Exception(msg)
+
+                    #print i, j, xi, eta, alpha, beta, vals[k], ref
+                    assert nanallclose(vals[k], ref, rtol=1e-12, atol=1e-12)
+
+    def test_linear_interpolation_outside_domain(self):
+        """Interpolation library sensibly handles values outside the domain
+        """
+
+        # Define pixel centers along each direction
+        x = [1.0, 2.0, 4.0]
+        y = [5.0, 9.0]
+
+        # Define ny by nx array with corresponding values
+        A = numpy.zeros((len(x), len(y)))
+
+        # Define values for each x, y pair as a linear function
+        for i in range(len(x)):
+            for j in range(len(y)):
+                A[i, j] = linear_function(x[i], y[j])
+
+        # Simple example first for debugging
+        xis = numpy.linspace(0.9, 4.0, 4)
+        etas = numpy.linspace(5, 9.1, 3)
+        points = axes2points(xis, etas)
+        refs = linear_function(points[:, 0], points[:, 1])
+
+        vals = interpolate2d(x, y, A, points, mode='linear',
+                             bounds_error=False)
+        msg = ('Length of interpolation points %i differs from length '
+               'of interpolated values %i' % (len(points), len(vals)))
+        assert len(points) == len(vals), msg
+        for i, (xi, eta) in enumerate(points):
+            if xi < x[0] or xi > x[-1] or eta < y[0] or eta > y[-1]:
+                assert numpy.isnan(vals[i])
+            else:
+                msg = ('Got %.15f for (%f, %f), expected %.15f'
+                       % (vals[i], xi, eta, refs[i]))
+                assert numpy.allclose(vals[i], refs[i],
+                                      rtol=1.0e-12, atol=1.0e-12), msg
+
+        # Try a range of combinations of points outside domain
+        # with error_bounds True
+        for lox in [x[0], x[0] - 1]:
+            for hix in [x[-1], x[-1] + 1]:
+                for loy in [y[0], y[0] - 1]:
+                    for hiy in [y[-1], y[-1] + 1]:
+
+                        # Then test that points outside domain can be handled
+                        xis = numpy.linspace(lox, hix, 4)
+                        etas = numpy.linspace(loy, hiy, 4)
+                        points = axes2points(xis, etas)
+
+                        if lox < x[0] or hix > x[-1] or \
+                                loy < x[0] or hiy > y[-1]:
+                            try:
+                                vals = interpolate2d(x, y, A, points,
+                                                     mode='linear',
+                                                     bounds_error=True)
+                            except Exception, e:
+                                pass
+                            else:
+                                msg = 'Should have raise bounds error'
+                                raise Exception(msg)
+
+        # Try a range of combinations of points outside domain with
+        # error_bounds False
+        for lox in [x[0], x[0] - 1, x[0] - 10]:
+            for hix in [x[-1], x[-1] + 1, x[-1] + 5]:
+                for loy in [y[0], y[0] - 1, y[0] - 10]:
+                    for hiy in [y[-1], y[-1] + 1, y[-1] + 10]:
+
+                        # Then test that points outside domain can be handled
+                        xis = numpy.linspace(lox, hix, 10)
+                        etas = numpy.linspace(loy, hiy, 10)
+                        points = axes2points(xis, etas)
+                        refs = linear_function(points[:, 0], points[:, 1])
+                        vals = interpolate2d(x, y, A, points,
+                                             mode='linear', bounds_error=False)
+
+                        assert len(points) == len(vals), msg
+                        for i, (xi, eta) in enumerate(points):
+                            if xi < x[0] or xi > x[-1] or\
+                                    eta < y[0] or eta > y[-1]:
+                                msg = 'Expected NaN for %f, %f' % (xi, eta)
+                                assert numpy.isnan(vals[i]), msg
+                            else:
+                                msg = ('Got %.15f for (%f, %f), expected '
+                                       '%.15f' % (vals[i], xi, eta, refs[i]))
+                                assert numpy.allclose(vals[i], refs[i],
+                                                      rtol=1.0e-12,
+                                                      atol=1.0e-12), msg
+
+
+if __name__ == '__main__':
+    suite = unittest.makeSuite(Test_interpolate, 'test')
+    runner = unittest.TextTestRunner(verbosity=2)
+    runner.run(suite)
