Repository: fatiando/verde
Pull request #46 (2018-05-20 08:42:32)

Biharmonic spline interpolation
Implement the biharmonic spline interpolation using Green's
functions of Sandwell (1987). The `Spline` gridder allows 
exact interpolation by fixing the forces on the data points or
approximate interpolation using damping or forces on a grid.
Approximate solutions allow the use of weights to take care
of outliers.


Patch:
diff --git a/.pylintrc b/.pylintrc
index a3e4f57..2064e80 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -50,7 +50,7 @@ confidence=
 # --enable=similarities". If you want to run only the classes checker, but have
 # no Warning level messages displayed, use"--disable=all --enable=classes
 # --disable=W"
-disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,import-star-module-level,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,useless-suppression,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init
+disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init,similarities
 
 # Enable the message, report, category or checker with the given id(s). You can
 # either give multiple identifier separated by comma (,) or put this option
@@ -164,10 +164,10 @@ ignore-comments=yes
 ignore-docstrings=yes
 
 # Ignore imports when computing similarities.
-ignore-imports=no
+ignore-imports=yes
 
 # Minimum lines number of a similarity.
-min-similarity-lines=4
+min-similarity-lines=3
 
 
 [LOGGING]
diff --git a/doc/index.rst b/doc/index.rst
index 65dd890..de6a3ef 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -135,3 +135,4 @@ terms of the :ref:`BSD 3-clause License <license>`.
     tutorials/index.rst
     api/index.rst
     license.rst
+    references.rst
diff --git a/doc/references.rst b/doc/references.rst
new file mode 100644
index 0000000..dacd1e7
--- /dev/null
+++ b/doc/references.rst
@@ -0,0 +1,4 @@
+References
+==========
+
+.. [Sandwell1987] Sandwell, D. T. (1987). Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data. Geophysical Research Letters, 14(2), 139–142. https://doi.org/10.1029/GL014i002p00139
diff --git a/examples/chain_trend.py b/examples/chain_trend.py
index 062300e..8e42cb3 100644
--- a/examples/chain_trend.py
+++ b/examples/chain_trend.py
@@ -36,7 +36,7 @@ coordinates = projection(lon, lat)
 # Create a chain that fits a 2nd degree trend to the anomaly data and then a
 # standard gridder to the residuals
 chain = vd.Chain([('trend', vd.Trend(degree=2)),
-                  ('gridder', vd.ScipyGridder())])
+                  ('spline', vd.Spline(damping=1e-8))])
 print("Chained estimator:", chain)
 # Calling 'fit' will automatically run the data through the chain
 chain.fit(coordinates, total_field)
@@ -47,7 +47,7 @@ grid_trend = chain.named_steps['trend'].grid()
 print("\nTrend grid:")
 print(grid_trend)
 
-grid_residual = chain.named_steps['gridder'].grid()
+grid_residual = chain.named_steps['spline'].grid()
 print("\nResidual grid:")
 print(grid_residual)
 
@@ -65,7 +65,7 @@ crs = ccrs.PlateCarree()
 
 plt.figure(figsize=(7, 5))
 ax = plt.axes(projection=ccrs.Mercator())
-ax.set_title("Chained trend and gridder")
+ax.set_title("Chained trend and spline")
 maxabs = np.max(np.abs([grid.total_field_anomaly.min(),
                         grid.total_field_anomaly.max()]))
 pc = ax.pcolormesh(grid.longitude, grid.latitude, grid.total_field_anomaly,
diff --git a/examples/spline.py b/examples/spline.py
new file mode 100644
index 0000000..f08da01
--- /dev/null
+++ b/examples/spline.py
@@ -0,0 +1,87 @@
+"""
+Dealing with outliers using splines
+===================================
+
+Biharmonic spline interpolation is based on estimating vertical forces acting
+on an elastic sheet that yield deformations in the sheet equal to the observed
+data.
+The results are equivalent to using ``verde.ScipyGridder(method='cubic')`` but
+the interpolation is usually slower.
+The advantage of using :class:`verde.Spline` is that you can assign weights to
+the data to deal with large uncertainties or outliers.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+# We need these two classes to set proper ticklabels for Cartopy maps
+from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
+import pyproj
+import numpy as np
+import verde as vd
+
+# We'll test this on the Baja California shipborne bathymetry data
+data = vd.datasets.fetch_baja_bathymetry()
+region = vd.get_region(data.longitude, data.latitude)
+
+# Before gridding, we need to decimate the data to avoid aliasing.
+spacing = 10/60
+lon, lat, bathymetry = vd.block_reduce(data.longitude, data.latitude,
+                                       data.bathymetry_m, reduction=np.median,
+                                       spacing=spacing)
+
+# Now we can add some outliers to test our spline
+outliers = np.random.RandomState(2).randint(0, bathymetry.size, size=20)
+bathymetry[outliers] += 5000
+print("Index of outliers:", outliers)
+
+# Project the data using pyproj so that we can use it as input for the gridder.
+# We'll set the latitude of true scale to the mean latitude of the data.
+projection = pyproj.Proj(proj='merc', lat_ts=data.latitude.mean())
+coordinates = projection(lon, lat)
+
+# Create an array of weights and set the weights for the outliers to a very low
+# value
+weights = np.ones_like(bathymetry)
+weights[outliers] = 1e-5
+
+# Now we can fit two splines to our data, one with weights and one without.
+# The default spline will fit the data exactly.
+spline = vd.Spline().fit(coordinates, bathymetry)
+# The weights only work if we use an approximate least-squares solution.
+# This means that the spline won't be exact on the data points.
+# The easiest way of doing this is to apply some damping regularization to
+# smooth the solution.
+spline_weights = vd.Spline(damping=1e-8).fit(coordinates, bathymetry, weights)
+
+# We'll make two geographic grids, one for each spline, to compare the results
+grids = [sp.grid(region=region, spacing=spacing, projection=projection,
+                 dims=['latitude', 'longitude'], data_names=['bathymetry'])
+         for sp in [spline, spline_weights]]
+
+# Now we can plot the two grids side by side on Mercator maps
+fig, axes = plt.subplots(1, 2, figsize=(9, 6),
+                         subplot_kw=dict(projection=ccrs.Mercator()))
+titles = ['No Weights', 'Weights']
+crs = ccrs.PlateCarree()
+for ax, title, grid in zip(axes, titles, grids):
+    ax.set_title(title)
+    pc = ax.pcolormesh(grid.longitude, grid.latitude, grid.bathymetry,
+                       transform=crs, vmax=0, vmin=data.bathymetry_m.min())
+    cb = plt.colorbar(pc, ax=ax, orientation='horizontal', pad=0.05)
+    cb.set_label('bathymetry [m]')
+    # Plot the land as a solid color
+    ax.add_feature(cfeature.LAND, edgecolor='black', zorder=2)
+    # Plot the locations of the decimated data
+    ax.plot(lon, lat, '.k', markersize=0.5, transform=crs)
+    # Plot the locations of the outliers
+    ax.plot(lon[outliers], lat[outliers], 'xk', transform=crs,
+            label='Outliers')
+    ax.set_extent(region, crs=crs)
+    # Set the proper ticks for a Cartopy map
+    ax.set_xticks(np.arange(-114, -105, 2), crs=crs)
+    ax.set_yticks(np.arange(21, 30, 2), crs=crs)
+    ax.xaxis.set_major_formatter(LongitudeFormatter())
+    ax.yaxis.set_major_formatter(LatitudeFormatter())
+ax.legend(loc='upper right')
+plt.tight_layout()
+plt.show()
diff --git a/verde/__init__.py b/verde/__init__.py
index 7ec54f6..ac19549 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -12,6 +12,7 @@ from .grid_math import block_reduce, distance_mask
 from .scipy_bridge import ScipyGridder
 from .trend import Trend, trend_jacobian
 from .chain import Chain
+from .spline import Spline, spline_jacobian
 
 
 # Get the version number through versioneer
diff --git a/verde/base.py b/verde/base.py
index f3cf26a..f935ab8 100644
--- a/verde/base.py
+++ b/verde/base.py
@@ -400,3 +400,39 @@ def check_data(data):
     if not isinstance(data, tuple):
         data = (data,)
     return data
+
+
+def check_fit_input(coordinates, data, weights):
+    """
+    Validate the inputs to the fit method of gridders.
+
+    Checks that the coordinates, data, and weights (if given) all have the same
+    shape.
+
+    Parameters
+    ----------
+    coordinates : tuple of arrays
+        Arrays with the coordinates of each data point. Should be in the
+        following order: (easting, northing, vertical, ...).
+    data : array
+        The data values of each data point.
+    weights : None or array
+        If not None, then the weights assigned to each data point.
+        Typically, this should be 1 over the data uncertainty squared.
+
+    Returns
+    -------
+    validated_inputs
+        The validated inputs in the same order. If weights are given, will
+        ravel the array before returning.
+
+    """
+    if any(coord.shape != data.shape for coord in coordinates):
+        raise ValueError(
+            "Coordinate and data arrays must have the same shape.")
+    if weights is not None:
+        if weights.shape != data.shape:
+            raise ValueError(
+                "Weights must have the same shape as the data array.")
+        weights = weights.ravel()
+    return coordinates, data, weights
diff --git a/verde/scipy_bridge.py b/verde/scipy_bridge.py
index 7dd2586..37d45ff 100644
--- a/verde/scipy_bridge.py
+++ b/verde/scipy_bridge.py
@@ -8,7 +8,7 @@ from sklearn.utils.validation import check_is_fitted
 from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator, \
     CloughTocher2DInterpolator
 
-from .base import BaseGridder
+from .base import BaseGridder, check_fit_input
 from .coordinates import get_region
 
 
@@ -33,13 +33,15 @@ class ScipyGridder(BaseGridder):
 
     Attributes
     ----------
+    interpolator_ : scipy interpolator class
+        An instance of the corresponding scipy interpolator class.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
         :meth:`~verde.ScipyGridder.grid` and
         :meth:`~verde.ScipyGridder.scatter` methods.
-    interpolator_ : scipy interpolator class
-        An instance of the corresponding scipy interpolator class.
+    residual_ : array
+        The difference between the input data and the interpolated values.
 
     """
 
@@ -93,6 +95,8 @@ class ScipyGridder(BaseGridder):
         if weights is not None:
             warn("{} does not support weights and they will be ignored."
                  .format(self.__class__.__name__))
+        coordinates, data, weights = check_fit_input(coordinates, data,
+                                                     weights)
         easting, northing = coordinates[:2]
         self.region_ = get_region(easting, northing)
         points = np.column_stack((np.ravel(easting), np.ravel(northing)))
diff --git a/verde/spline.py b/verde/spline.py
index e9d9a4e..270e0c0 100644
--- a/verde/spline.py
+++ b/verde/spline.py
@@ -1,14 +1,227 @@
 """
-Spline interpolation using Green's functions.
-
-Estimates the intensity of vertical forces acting on an elastic line or sheet.
-The interpolated data are the deflections of this sheet by the forces.
+Biharmonic splines in 2D.
 """
+from warnings import warn
+
 import numpy as np
+from sklearn.utils.validation import check_is_fitted
+from sklearn.preprocessing import StandardScaler
+from sklearn.linear_model import LinearRegression, Ridge
+
+from .base import BaseGridder, check_fit_input
+from .coordinates import grid_coordinates, get_region
+
+
+class Spline(BaseGridder):
+    r"""
+    Biharmonic spline interpolation using Green's functions.
+
+    Implements the 2D splines of [Sandwell1987]_. The Green's function for the
+    spline corresponds to the elastic deflection of a thin sheet subject to a
+    vertical force. For an observation point at the origin and a force at the
+    coordinates given by the vector :math:`\mathbf{x}`, the Green's function
+    is:
+
+    .. math::
+
+        g(\mathbf{x}) = \|\mathbf{x}\|^2 \left(\log \|\mathbf{x}\| - 1\right)
+
+    In practice, this function is not defined for data points that coincide
+    with a force. To prevent this, a fudge factor is added to
+    :math:`\|\mathbf{x}\|`.
+
+    The interpolation is performed by estimating forces that produce
+    deflections that fit the observed data (using least-squares). Then, the
+    interpolated points can be evaluated at any location.
+
+    By default, the forces will be placed at the same points as the input data
+    given to :meth:`~verde.Spline.fit`. This configuration provides an exact
+    solution on top of the data points. However, this solution can be unstable
+    for certain configurations of data points.
+
+    Approximate (and more stable) solutions can be obtained by applying damping
+    regularization to smooth the estimated forces (and interpolated values) or
+    by distributing the forces on a regular grid instead of at the data points.
+
+    Data weights can be used during fitting but only have an any effect when
+    using the approximate solutions.
+
+    The Jacobian (design, sensitivity, feature, etc) matrix for the spline
+    is normalized using :class:`sklearn.preprocessing.StandardScaler` without
+    centering the mean so that the transformation can be undone in the
+    estimated forces.
+
+    Parameters
+    ----------
+    fudge : float
+        The positive fudge factor applied to the Green's function to avoid
+        singularities.
+    damping : None or float
+        The positive damping regularization parameter. Controls how much
+        smoothness is imposed on the estimated forces. If None, no
+        regularization is used.
+    shape : None or tuple
+        If not None, then should be the shape of the regular grid of forces.
+        See :func:`verde.grid_coordinates` for details.
+    spacing : None or float or tuple
+        If not None, then should be the spacing of the regular grid of forces.
+        See :func:`verde.grid_coordinates` for details.
+
+    Attributes
+    ----------
+    forces_ : array
+        The estimated forces that fit the observed data.
+    force_coords_ : tuple of arrays
+        The easting and northing coordinates of the forces.
+    region_ : tuple
+        The boundaries (``[W, E, S, N]``) of the data used to fit the
+        interpolator. Used as the default region for the
+        :meth:`~verde.Spline.grid` and :meth:`~verde.Spline.scatter` methods.
+    residual_ : array
+        The difference between the input data and the predicted spline values.
+
+    See also
+    --------
+    verde.spline_jacobian : Make the Jacobian matrix for the biharmonic spline
+
+    """
+
+    def __init__(self, fudge=1e-5, damping=None, shape=None, spacing=None):
+        self.damping = damping
+        self.shape = shape
+        self.spacing = spacing
+        self.fudge = fudge
+
+    def fit(self, coordinates, data, weights=None):
+        """
+        Fit the biharmonic spline to the given data.
+
+        The data region is captured and used as default for the
+        :meth:`~verde.Spline.grid` and :meth:`~verde.Spline.scatter` methods.
+
+        All input arrays must have the same shape.
+
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting
+            and northing will be used, all subsequent coordinates will be
+            ignored.
+        data : array
+            The data values of each data point.
+        weights : None or array
+            If not None, then the weights assigned to each data point.
+            Typically, this should be 1 over the data uncertainty squared.
+
+        Returns
+        -------
+        self
+            Returns this estimator instance for chaining operations.
+
+        """
+        coordinates, data, weights = check_fit_input(coordinates, data,
+                                                     weights)
+        exact_solution = all(i is None
+                             for i in [self.damping, self.spacing, self.shape])
+        if weights is not None and exact_solution:
+            warn(' '.join([
+                "Weights are ignored for the exact spline solution.",
+                "Use damping or specify a shape/spacing for the forces."]))
+        easting, northing = coordinates[:2]
+        self.region_ = get_region(easting, northing)
+        # Set the force positions. If no shape or spacing are given, then they
+        # will coincide with the data points.
+        if self.shape is None and self.spacing is None:
+            self.force_coords_ = (easting.copy(), northing.copy())
+        else:
+            self.force_coords_ = grid_coordinates(
+                self.region_, shape=self.shape, spacing=self.spacing)
+        jac = spline_jacobian(easting, northing, self.force_coords_[0],
+                              self.force_coords_[1], self.fudge)
+        scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
+        jac = scaler.fit_transform(jac)
+        if self.damping is None:
+            regr = LinearRegression(fit_intercept=False, normalize=False)
+        else:
+            regr = Ridge(alpha=self.damping, fit_intercept=False,
+                         normalize=False)
+        regr.fit(jac, data.ravel(), sample_weight=weights)
+        self.residual_ = data - jac.dot(regr.coef_).reshape(data.shape)
+        self.force_ = regr.coef_/scaler.scale_
+        return self
+
+    def predict(self, coordinates):
+        """
+        Evaluate the estimated spline on the given set of points.
 
+        Requires a fitted estimator (see :meth:`~verde.Trend.fit`).
 
-def spline1d_green():
+        Parameters
+        ----------
+        coordinates : tuple of arrays
+            Arrays with the coordinates of each data point. Should be in the
+            following order: (easting, northing, vertical, ...). Only easting
+            and northing will be used, all subsequent coordinates will be
+            ignored.
+
+        Returns
+        -------
+        data : array
+            The data values evaluated on the given points.
+
+        """
+        easting, northing = coordinates[:2]
+        check_is_fitted(self, ['force_', 'force_coords_'])
+        jac = spline_jacobian(easting.ravel(), northing.ravel(),
+                              self.force_coords_[0], self.force_coords_[1],
+                              self.fudge)
+        shape = np.broadcast(easting, northing).shape
+        return jac.dot(self.force_).reshape(shape)
+
+
+def spline_jacobian(easting, northing, force_easting, force_northing,
+                    fudge=1e-5):
     """
-    The Green's function for the 1D spline.
+    Make the Jacobian matrix for the 2D biharmonic spline.
+
+    Each column of the Jacobian is the Green's function for a single force
+    evaluated on all observation points.
+
+    Parameters
+    ----------
+    easting : array
+        The values of the West-East coordinates of each data point.
+    northing : array
+        The values of the South-North coordinates of each data point.
+    force_easting : array
+        The values of the West-East coordinates of each force.
+    force_northing : array
+        The values of the South-North coordinates of each force.
+    fudge : float
+        The positive fudge factor applied to the Green's function to avoid
+        singularities.
+
+    Returns
+    -------
+    jacobian : 2D array
+        The (n_data, n_forces) Jacobian matrix.
+
+    See also
+    --------
+    verde.Spline : Biharmonic spline gridder
+
     """
-    return np.ones(10)
+    easting = np.atleast_1d(easting)
+    northing = np.atleast_1d(northing)
+    force_easting = np.atleast_1d(force_easting)
+    force_northing = np.atleast_1d(force_northing)
+    # Reshaping the data to a column vector will automatically build a
+    # distance matrix between each data point and force.
+    distance = np.hypot(
+        easting.reshape((easting.size, 1)) - force_easting.ravel(),
+        northing.reshape((northing.size, 1)) - force_northing.ravel())
+    # The fudge factor helps avoid singular matrices when the force and
+    # computation point are too close
+    distance += fudge
+    return (distance**2)*(np.log(distance) - 1)
diff --git a/verde/trend.py b/verde/trend.py
index e589dc6..aacf16d 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -6,7 +6,7 @@ from sklearn.utils.validation import check_is_fitted
 from sklearn.preprocessing import StandardScaler
 from sklearn.linear_model import LinearRegression
 
-from .base import BaseGridder
+from .base import BaseGridder, check_fit_input
 from .coordinates import get_region
 
 
@@ -26,6 +26,18 @@ class Trend(BaseGridder):
     degree : int
         The degree of the polynomial. Must be >= 1.
 
+    Attributes
+    ----------
+    coef_ : array
+        The estimated polynomial coefficients that fit the observed data.
+    region_ : tuple
+        The boundaries (``[W, E, S, N]``) of the data used to fit the
+        interpolator. Used as the default region for the
+        :meth:`~verde.Trend.grid` and :meth:`~verde.Trend.scatter` methods.
+    residual_ : array
+        The difference between the input data and the predicted polynomial
+        trend.
+
     Examples
     --------
 
@@ -51,6 +63,10 @@ class Trend(BaseGridder):
     >>> print('{:.2f}'.format(trend_out.residual_[2, 2]))
     500.00
 
+    See also
+    --------
+    verde.trend_jacobian ; Make the Jacobian matrix for 2D polynomial
+
     """
 
     def __init__(self, degree):
@@ -80,19 +96,13 @@ class Trend(BaseGridder):
 
         Returns
         -------
-        self : verde.Trend
+        self
             Returns this estimator instance for chaining operations.
 
         """
+        coordinates, data, weights = check_fit_input(coordinates, data,
+                                                     weights)
         easting, northing = coordinates[:2]
-        if easting.shape != northing.shape or easting.shape != data.shape:
-            raise ValueError(
-                "Coordinate and data arrays must have the same shape.")
-        if weights is not None:
-            if weights.shape != data.shape:
-                raise ValueError(
-                    "Weights must have the same shape as the data array.")
-            weights = weights.ravel()
         self.region_ = get_region(easting, northing)
         jac = trend_jacobian(easting, northing, degree=self.degree,
                              dtype=data.dtype)
@@ -189,7 +199,7 @@ def trend_jacobian(easting, northing, degree, dtype='float64'):
     Returns
     -------
     jacobian : 2D array
-        The Jacobian matrix.
+        The (n_data, n_coefficients) Jacobian matrix.
 
     Examples
     --------
@@ -210,6 +220,10 @@ def trend_jacobian(easting, northing, degree, dtype='float64'):
      [ 1  3 -2  9 -6  4]
      [ 1  4 -1 16 -4  1]]
 
+    See also
+    --------
+    verde.Trend ; Polynomial trend estimator
+
     """
     if easting.shape != northing.shape:
         raise ValueError("Coordinate arrays must have the same shape.")


Test patch:
diff --git a/verde/tests/test_base.py b/verde/tests/test_base.py
index 02f1eb2..ab4953e 100644
--- a/verde/tests/test_base.py
+++ b/verde/tests/test_base.py
@@ -6,7 +6,8 @@ import pytest
 import numpy as np
 import numpy.testing as npt
 
-from ..base import get_dims, get_data_names, get_region, BaseGridder
+from ..base import get_dims, get_data_names, get_region, BaseGridder, \
+    check_fit_input
 from ..coordinates import grid_coordinates, scatter_points
 
 
@@ -182,3 +183,20 @@ def test_basegridder_projection():
     npt.assert_allclose(scat.scalars, data)
     npt.assert_allclose(prof.scalars,
                         angular*coordinates_true[0][0, :] + linear)
+
+
+def test_check_fit_input_fails_coordinates():
+    "Test the failing conditions for check_fit_input"
+    coords = (np.arange(20), np.arange(20))
+    data = np.arange(30)
+    with pytest.raises(ValueError):
+        check_fit_input(coords, data, weights=None)
+
+
+def test_check_fit_input_fails_weights():
+    "Test the failing conditions for check_fit_input"
+    data = np.arange(20)
+    coords = (data, data)
+    weights = np.arange(30)
+    with pytest.raises(ValueError):
+        check_fit_input(coords, data, weights)
diff --git a/verde/tests/test_spline.py b/verde/tests/test_spline.py
index ce3f25d..267b8d4 100644
--- a/verde/tests/test_spline.py
+++ b/verde/tests/test_spline.py
@@ -1,12 +1,93 @@
 """
-Test the spline interpolation.
+Test the biharmonic splines
 """
+import warnings
+
 import numpy as np
 import numpy.testing as npt
 
-from ..spline import spline1d_green
+from ..spline import Spline
+from ..datasets.synthetic import CheckerBoard
+
+
+def test_spline():
+    "See if the exact spline solution."
+    region = (100, 500, -800, -700)
+    synth = CheckerBoard().fit(region=region)
+    data = synth.scatter(size=1500, random_state=1)
+    coords = (data.easting, data.northing)
+    # The interpolation should be perfect on top of the data points
+    spline = Spline().fit(coords, data.scalars)
+    npt.assert_allclose(spline.predict(coords), data.scalars, rtol=1e-5)
+    npt.assert_allclose(spline.residual_, 0, atol=1e-5)
+    # There should be 1 force per data point
+    assert data.scalars.size == spline.force_.size
+    npt.assert_allclose(spline.force_coords_, coords)
+    shape = (5, 5)
+    region = (270, 320, -770, -720)
+    # Tolerance needs to be kind of high to allow for error due to small
+    # dataset
+    npt.assert_allclose(spline.grid(region, shape=shape).scalars,
+                        synth.grid(region, shape=shape).scalars,
+                        rtol=5e-2)
+
+
+def test_spline_region():
+    "See if the region is gotten from the data is correct."
+    region = (1000, 5000, -8000, -6000)
+    grid = CheckerBoard().fit(region=region).grid(shape=(10, 10))
+    coords = np.meshgrid(grid.easting, grid.northing)
+    grd = Spline().fit(coords, grid.scalars.values)
+    npt.assert_allclose(grd.region_, region)
+
+
+def test_spline_damping():
+    "Test the spline solution with a bit of damping"
+    region = (100, 500, -800, -700)
+    synth = CheckerBoard().fit(region=region)
+    data = synth.scatter(size=1200, random_state=1)
+    coords = (data.easting, data.northing)
+    # The interpolation should be close on top of the data points
+    spline = Spline(damping=1e-15).fit(coords, data.scalars)
+    npt.assert_allclose(spline.residual_.mean(), 0, atol=1e-2)
+    npt.assert_allclose(spline.residual_.std(), 0, atol=1e-2)
+    shape = (5, 5)
+    region = (270, 320, -770, -720)
+    # Tolerance needs to be kind of high to allow for error due to small
+    # dataset
+    npt.assert_allclose(spline.grid(region, shape=shape).scalars,
+                        synth.grid(region, shape=shape).scalars,
+                        rtol=5e-2)
+
+
+def test_spline_grid():
+    "Test the spline solution with forces on a grid"
+    region = (1000, 5000, -8000, -7000)
+    synth = CheckerBoard().fit(region=region)
+    data = synth.scatter(size=1500, random_state=1)
+    coords = (data.easting, data.northing)
+    # The interpolation should be close on top of the data points
+    spline = Spline(shape=(35, 35)).fit(coords, data.scalars)
+    npt.assert_allclose(spline.residual_.mean(), 0, atol=1e-2)
+    npt.assert_allclose(spline.residual_.std(), 0, atol=2)
+    assert spline.force_.size == 35**2
+    shape = (3, 3)
+    region = (2700, 3200, -7700, -7200)
+    # Tolerance needs to be kind of high to allow for error due to small
+    # dataset
+    npt.assert_allclose(spline.grid(region, shape=shape).scalars,
+                        synth.grid(region, shape=shape).scalars,
+                        rtol=5e-2)
 
 
-def test_spline1d_green():
-    "Test the spline 1D Greens function against known values"
-    npt.assert_allclose(spline1d_green(), np.ones(10))
+def test_spline_warns():
+    "Check that a warning is issued when using weights and the exact solution."
+    data = CheckerBoard().fit().scatter(random_state=100)
+    weights = np.ones_like(data.scalars)
+    grd = Spline()
+    msg = "Weights are ignored for the exact spline solution"
+    with warnings.catch_warnings(record=True) as warn:
+        grd.fit((data.easting, data.northing), data.scalars, weights=weights)
+        assert len(warn) == 1
+        assert issubclass(warn[-1].category, UserWarning)
+        assert str(warn[-1].message).split('.')[0] == msg
diff --git a/verde/tests/test_trend.py b/verde/tests/test_trend.py
index 70ac122..fa9a1f1 100644
--- a/verde/tests/test_trend.py
+++ b/verde/tests/test_trend.py
@@ -42,23 +42,6 @@ def test_trend_weights(simple_model):
     npt.assert_allclose(trend.predict(coords), data)
 
 
-def test_trend_fails_coordinates():
-    "Test the failing conditions for the trend estimator"
-    coords = (np.arange(20), np.arange(20))
-    data = np.arange(30)
-    with pytest.raises(ValueError):
-        Trend(degree=10).fit(coords, data)
-
-
-def test_trend_fails_weights():
-    "Test the failing conditions for the trend estimator"
-    data = np.arange(20)
-    coords = (data, data)
-    weights = np.arange(30)
-    with pytest.raises(ValueError):
-        Trend(degree=10).fit(coords, data, weights)
-
-
 def test_polynomial_combinations_fails():
     "Test failing conditions for the combinations"
     with pytest.raises(ValueError):