Repository: fatiando/verde
Pull request #67 (2018-05-27 23:39:33)

Implement the Vector2D gridder
Implements the method of Sandwell and Wessel (2016) for gridding
2-component vector data using elastic Green's functions in 
`verde.Vector2D`.
Extends the method by using ridge regression instead of SVD
and adding the option of laying the forces out in a grid (which is 
required for KFold cross-validation because otherwise the model
would change for different data folds).

Refactors `verde.Spline` so that it can be used as a base class for
`verde.Vector2D` to avoid code repetition.


leouieda: TODO:

- [ ] Write an example using the simple classic solution without regularization
- [ ] Split out the cross-validation stuff into a separate PR
- [ ] Write tests for Vector2D
- [ ] Write the docstrings for Vector2D

Patch:
diff --git a/doc/install.rst b/doc/install.rst
index 73ea813..80d21c5 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -25,10 +25,13 @@ Dependencies
 * `scikit-learn <http://scikit-learn.org/>`__
 * `requests <http://docs.python-requests.org/>`__
 
-Most of the examples in the `Gallery </gallery/index.html>`__ and
-`Tutorials </tutorials/index.html>`__
-also use `matplotlib <https://matplotlib.org/>`__ and
-`cartopy <https://scitools.org.uk/cartopy/>`__ for plotting.
+Most of the examples in the :ref:`gallery` and :ref:`tutorials` also use:
+
+* `matplotlib <https://matplotlib.org/>`__
+* `cartopy <https://scitools.org.uk/cartopy/>`__ for plotting maps
+* `pyproj <https://jswhit.github.io/pyproj/>`__ for cartographic projections
+* `dask <https://dask.pydata.org/>`__ for parallelism
+
 
 
 Installing the latest development version
@@ -50,10 +53,7 @@ Testing your install
 
 We ship a full test suite with the package.
 You can run our tests after you install it but you will need to install `pytest
-<https://docs.pytest.org/>`__ first::
-
-    conda install -c conda-forge pytest
-
+<https://docs.pytest.org/>`__ first.
 After that, you can test your installation by running the following inside a
 Python interpreter::
 
diff --git a/doc/references.rst b/doc/references.rst
index dacd1e7..38a4bbc 100644
--- a/doc/references.rst
+++ b/doc/references.rst
@@ -2,3 +2,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
+.. [SandwellWessel2016] Sandwell, D. T., & Wessel, P. (2016). Interpolation of 2-D vector data using constraints from elasticity. Geophysical Research Letters, 43(20), 2016GL070340. https://doi.org/10.1002/2016GL070340
diff --git a/environment.yml b/environment.yml
index 7c9f23d..26e8e7e 100644
--- a/environment.yml
+++ b/environment.yml
@@ -12,6 +12,8 @@ dependencies:
     - scikit-learn
     - requests
     # Development requirements
+    - dask
+    - pyproj
     - matplotlib
     - cartopy
     - jupyter
diff --git a/examples/README.txt b/examples/README.txt
index 0c4b5cc..c1a8270 100644
--- a/examples/README.txt
+++ b/examples/README.txt
@@ -1,3 +1,5 @@
+.. _gallery:
+
 Gallery
 =======
 
diff --git a/examples/vector2d.py b/examples/vector2d.py
new file mode 100644
index 0000000..8d862f3
--- /dev/null
+++ b/examples/vector2d.py
@@ -0,0 +1,113 @@
+"""
+Gridding 2-component vector data
+================================
+
+One way of gridding vector data would be grid each component separately using
+:class:`verde.Spline`. Alternatively, :class:`verde.Vector2D` can grid two
+components simultaneously in a way that couples them through elastic
+deformation theory. This is particularly suited, though not exclusive, to data
+that represent elastic/semi-elastic deformation, like horizontal GPS
+velocities.
+"""
+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 numpy as np
+import pyproj
+import verde as vd
+
+
+# Fetch the GPS data from the U.S. West coast. We'll grid only the horizontal
+# components of the velocities
+data = vd.datasets.fetch_california_gps()
+coordinates = (data.longitude.values, data.latitude.values)
+region = vd.get_region(coordinates)
+# Use a Mercator projection because Vector2D is a Cartesian gridder
+projection = pyproj.Proj(proj='merc', lat_ts=data.latitude.mean())
+
+# Split the data into a training and testing set. We'll fit the gridder on the
+# training set and use the testing set to evaluate how well the gridder is
+# performing.
+train, test = vd.train_test_split(projection(*coordinates),
+                                  (data.velocity_east, data.velocity_north),
+                                  random_state=0)
+
+# We'll make a 20 arc-minute grid in the end.
+spacing = 20/60
+
+# Chain together a blocked mean to avoid aliasing, a polynomial trend to take
+# care of the increase toward the coast, and finally the vector gridder using
+# Poisson's ratio 0.5 to couple the two horizontal components.
+chain = vd.Chain([('mean', vd.BlockReduce(np.mean, spacing*111e3)),
+                  ('trend', vd.VectorTrend(degree=5)),
+                  ('spline', vd.Vector2D(poisson=0.5))])
+# Fit on the training data
+chain.fit(*train)
+# And score on the testing data. The best possible score is 1, meaning a
+# perfect prediction of the test data.
+score = chain.score(*test)
+print("Cross-validation R^2 score: {:.2f}".format(score))
+
+# Interpolate our horizontal GPS velocities onto a regular geographic grid
+grid = chain.grid(region, spacing=spacing, projection=projection,
+                  dims=['latitude', 'longitude'])
+# Mask the data that are far from the observation points
+mask = vd.distance_mask(np.meshgrid(grid.longitude, grid.latitude),
+                        (data.longitude, data.latitude),
+                        maxdist=0.5)
+grid = grid.where(~mask)
+
+# Calculate residuals between the predictions and the original input data.
+# Even though we aren't using regularization or regularly distributed forces,
+# the prediction won't be perfect because of the BlockReduce operation. We fit
+# the gridder on the reduced observations, not the original data.
+predicted = chain.predict(projection(*coordinates))
+residuals = (data.velocity_east - predicted[0],
+             data.velocity_north - predicted[1])
+
+# Make maps of the original velocities, the gridded velocities, and the
+# residuals
+
+
+def setup_map(ax, title):
+    "Setup a map with land and ocean polygons and proper tick labels"
+    ax.set_title(title)
+    # Setup the map ticks
+    ax.set_xticks(np.arange(-124, -115, 2), crs=crs)
+    ax.set_yticks(np.arange(33, 42, 2), crs=crs)
+    ax.xaxis.set_major_formatter(LongitudeFormatter())
+    ax.yaxis.set_major_formatter(LatitudeFormatter())
+    ax.add_feature(cfeature.LAND)
+    ax.add_feature(cfeature.OCEAN)
+    # Extend the plot region to avoid cutting off some vectors at the borders
+    ax.set_extent(vd.pad_region(region, pad=0.5), crs=crs)
+
+
+fig, axes = plt.subplots(1, 2, figsize=(12, 8),
+                         subplot_kw=dict(projection=ccrs.Mercator()))
+crs = ccrs.PlateCarree()
+# Plot the observed data and the residuals
+ax = axes[0]
+tmp = ax.quiver(data.longitude.values, data.latitude.values,
+                data.velocity_east.values, data.velocity_north.values,
+                scale=0.3, transform=crs, width=0.001, label="Velocities")
+ax.quiverkey(tmp, 0.13, 0.18, 0.05, label='0.05 m/yr', coordinates='figure')
+ax.quiver(data.longitude.values, data.latitude.values,
+          residuals[0].values, residuals[1].values,
+          scale=0.3, transform=crs, color='r', width=0.001, label="Residuals")
+setup_map(ax, 'GPS horizontal velocities')
+ax.legend(loc="lower left")
+# Plot the gridded data and the residuals
+ax = axes[1]
+tmp = ax.quiver(grid.longitude.values, grid.latitude.values,
+                grid.east_component.values, grid.north_component.values,
+                scale=0.3, transform=crs, width=0.002)
+ax.quiverkey(tmp, 0.63, 0.18, 0.05, label='0.05 m/yr', coordinates='figure')
+ax.quiver(data.longitude.values, data.latitude.values,
+          residuals[0].values, residuals[1].values,
+          scale=0.3, transform=crs, color='r', width=0.001)
+setup_map(ax, 'Gridded velocities')
+plt.tight_layout()
+plt.show()
diff --git a/requirements.txt b/requirements.txt
index e60f2f8..cb80f6b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,6 +5,8 @@ xarray
 scikit-learn
 requests
 # Development requirements
+dask
+pyproj
 matplotlib
 cartopy
 jupyter
diff --git a/tutorials/README.txt b/tutorials/README.txt
index 0c7e28c..4802055 100644
--- a/tutorials/README.txt
+++ b/tutorials/README.txt
@@ -1,2 +1,4 @@
+.. _tutorials:
+
 Tutorials
 =========
diff --git a/verde/__init__.py b/verde/__init__.py
index 9b1fa65..0aaa593 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -16,6 +16,7 @@ from .trend import Trend, trend_jacobian, VectorTrend
 from .chain import Chain
 from .spline import Spline, spline_jacobian
 from .model_selection import cross_val_score, train_test_split
+from .vector import Vector2D, vector2d_jacobian
 
 
 # Get the version number through versioneer
diff --git a/verde/base.py b/verde/base.py
index 4907b68..7640e54 100644
--- a/verde/base.py
+++ b/verde/base.py
@@ -214,7 +214,7 @@ class BaseGridder(BaseEstimator):
 
         """
         coordinates, data, weights = check_fit_input(coordinates, data,
-                                                     weights, ravel=False)
+                                                     weights, unpack=False)
         pred = check_data(self.predict(coordinates))
         result = np.mean(
             [r2_score(datai, predi, sample_weight=weighti)
@@ -529,12 +529,12 @@ def check_data(data):
     return data
 
 
-def check_fit_input(coordinates, data, weights, ravel=True):
+def check_fit_input(coordinates, data, weights, unpack=True):
     """
     Validate the inputs to the fit method of gridders.
 
     Checks that the coordinates, data, and weights (if given) all have the same
-    shape.
+    shape. Weights arrays are raveled.
 
     Parameters
     ----------
@@ -549,10 +549,10 @@ def check_fit_input(coordinates, data, weights, ravel=True):
         Typically, this should be 1 over the data uncertainty squared.
         If the data has multiple components, the weights have the same number
         of components.
-    ravel : bool
+    unoack : bool
         If False, data and weights will be tuples always. If they are single
         arrays, then they will be returned as a 1-element tuple. If True, will
-        only return tuples if there are more than 1 array in each.
+        unpack the tuples if there is only 1 array in each.
 
     Returns
     -------
@@ -577,7 +577,7 @@ def check_fit_input(coordinates, data, weights, ravel=True):
         weights = tuple(i.ravel() for i in weights)
     else:
         weights = tuple([None]*len(data))
-    if ravel:
+    if unpack:
         if len(weights) == 1:
             weights = weights[0]
         if len(data) == 1:
diff --git a/verde/blockreduce.py b/verde/blockreduce.py
index 789dfa9..2c7a641 100644
--- a/verde/blockreduce.py
+++ b/verde/blockreduce.py
@@ -199,7 +199,7 @@ class BlockReduce(BaseEstimator):
 
         """
         coordinates, data, weights = check_fit_input(coordinates, data,
-                                                     weights, ravel=False)
+                                                     weights, unpack=False)
         blocks, labels = block_split(coordinates, self.spacing, self.adjust,
                                      self.region)
         if any(w is None for w in weights):
@@ -390,7 +390,7 @@ class BlockMean(BlockReduce):
 
         """
         coordinates, data, weights = check_fit_input(coordinates, data,
-                                                     weights, ravel=False)
+                                                     weights, unpack=False)
         if any(w is None for w in weights) and self.uncertainty:
             raise ValueError(
                 "Weights are required for uncertainty propagation."
diff --git a/verde/chain.py b/verde/chain.py
index 73c66ce..7e02f4b 100644
--- a/verde/chain.py
+++ b/verde/chain.py
@@ -121,5 +121,5 @@ class Chain(BaseGridder):
                 for i, pred in enumerate(predicted):
                     result[i] += pred
         if len(result) == 1:
-            result = result[0]
-        return result
+            return result[0]
+        return tuple(result)
diff --git a/verde/model_selection.py b/verde/model_selection.py
index fd8eda4..d167a03 100644
--- a/verde/model_selection.py
+++ b/verde/model_selection.py
@@ -95,7 +95,7 @@ def train_test_split(coordinates, data, weights=None, **kwargs):
       (None,)
 
     """
-    args = check_fit_input(coordinates, data, weights, ravel=False)
+    args = check_fit_input(coordinates, data, weights, unpack=False)
     ndata = args[1][0].size
     indices = np.arange(ndata)
     split = next(ShuffleSplit(n_splits=1, **kwargs).split(indices))
@@ -178,7 +178,7 @@ def cross_val_score(estimator, coordinates, data, weights=None, cv=None,
 
     """
     coordinates, data, weights = check_fit_input(coordinates, data, weights,
-                                                 ravel=False)
+                                                 unpack=False)
     if client is None:
         client = DummyClient()
     if cv is None:
diff --git a/verde/spline.py b/verde/spline.py
index 6706065..c137f39 100644
--- a/verde/spline.py
+++ b/verde/spline.py
@@ -127,43 +127,75 @@ class Spline(BaseGridder):
         """
         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."]))
+        self._check_weighted_exact_solution(weights)
+        # Capture the data region to use as a default when gridding.
         self.region_ = get_region(coordinates[:2])
-        # 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_ = tuple(i.copy() for i in coordinates[:2])
-        else:
-            # Set the force grid region to the data region if none is given.
-            # This will only happen once so the model won't change during later
-            # fits.
-            if self.region is None:
-                self.region = self.region_
-            self.force_coords_ = grid_coordinates(self.region,
-                                                  shape=self.shape,
-                                                  spacing=self.spacing)
-        jac = spline_jacobian(coordinates[:2], self.force_coords_, self.fudge)
+        self.force_coords_ = self._get_force_coordinates(coordinates)
+        jacobian = spline_jacobian(coordinates[:2], self.force_coords_,
+                                   self.fudge)
+        self.force_ = self._estimate_forces(jacobian, data, weights)
+        return self
+
+    def _estimate_forces(self, jacobian, data, weights):
+        """
+        Estimate forces that fit the data using least-squares. Scales the
+        Jacobian matrix to have unit standard deviation. This helps balance the
+        regularization and the difference between forces.
+        """
+        if jacobian.shape[0] < jacobian.shape[1]:
+            warn(" ".join([
+                "Under-determined problem detected",
+                "(ndata, nparams)={}.".format(jacobian.shape),
+                "Configuration of forces:",
+                "spacing={} shape={}".format(self.spacing, self.shape)]))
         scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
-        jac = scaler.fit_transform(jac)
+        jacobian = scaler.fit_transform(jacobian)
         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.force_ = regr.coef_/scaler.scale_
-        return self
+        regr.fit(jacobian, data.ravel(), sample_weight=weights)
+        # Undo the scaling so that we can use forces on the unscaled Jacobian
+        # later on.
+        forces = regr.coef_/scaler.scale_
+        return forces
+
+    def _get_force_coordinates(self, coordinates):
+        """
+        Generate arrays for the force coordinates depending on the
+        configuration of the spline.
+        """
+        # Set the force positions. If no shape or spacing are given, then they
+        # will coincide with the data points. This will only happen once so the
+        # model won't change during later fits.
+        if self.shape is None and self.spacing is None:
+            coords = tuple(i.copy() for i in coordinates[:2])
+        else:
+            if self.region is None:
+                self.region = get_region(coordinates)
+            coords = grid_coordinates(self.region, shape=self.shape,
+                                      spacing=self.spacing)
+        return coords
+
+    def _check_weighted_exact_solution(self, weights):
+        """
+        If a weighted exact solution is requested, warn the user that it won't
+        work.
+        """
+        # Check if we're using weights with the exact solution and warn the
+        # user that it won't have any effect.
+        exact = all(i is None for i in [self.damping, self.spacing,
+                                        self.shape])
+        if weights is not None and exact:
+            warn("Weights are ignored for the exact solution. "
+                 "Use damping or specify a shape/spacing for the forces.")
 
     def predict(self, coordinates):
         """
         Evaluate the estimated spline on the given set of points.
 
-        Requires a fitted estimator (see :meth:`~verde.Trend.fit`).
+        Requires a fitted estimator (see :meth:`~verde.Spline.fit`).
 
         Parameters
         ----------
@@ -189,6 +221,8 @@ def spline_jacobian(coordinates, force_coordinates, fudge=1e-5):
     """
     Make the Jacobian matrix for the 2D biharmonic spline.
 
+    Follows [Sandwell1987]_.
+
     Each column of the Jacobian is the Green's function for a single force
     evaluated on all observation points.
 
@@ -216,12 +250,9 @@ def spline_jacobian(coordinates, force_coordinates, fudge=1e-5):
     verde.Spline : Biharmonic spline gridder
 
     """
-    force_easting, force_northing = force_coordinates[:2]
-    easting, northing = coordinates[:2]
-    easting = np.atleast_1d(easting).ravel()
-    northing = np.atleast_1d(northing).ravel()
-    force_easting = np.atleast_1d(force_easting)
-    force_northing = np.atleast_1d(force_northing)
+    force_easting, force_northing = (np.atleast_1d(i).ravel()
+                                     for i in force_coordinates[:2])
+    easting, northing = (np.atleast_1d(i).ravel() for i in coordinates[:2])
     # Reshaping the data to a column vector will automatically build a
     # distance matrix between each data point and force.
     distance = np.hypot(
diff --git a/verde/trend.py b/verde/trend.py
index bf5c6b1..35f8de9 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -252,7 +252,7 @@ class VectorTrend(BaseGridder):
 
     Attributes
     ----------
-    components_ : list
+    component_ : list
         List of the fitted :class:`verde.Trend` on each component of the data.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
diff --git a/verde/vector.py b/verde/vector.py
new file mode 100644
index 0000000..d77c88b
--- /dev/null
+++ b/verde/vector.py
@@ -0,0 +1,244 @@
+"""
+Vector gridding using elasticity Green's functions from Sandwell and Wessel
+(2016).
+"""
+import numpy as np
+from sklearn.utils.validation import check_is_fitted
+
+from .base import check_fit_input
+from .spline import Spline
+from .coordinates import get_region
+
+
+class Vector2D(Spline):
+    r"""
+    Elastically coupled interpolation of 2-component vector data.
+
+    Uses the Green's functions based on elastic deformation from
+    [SandwellWessel2016]_. The interpolation is done by estimating point forces
+    that generate an elastic deformation that fits the observed vector data.
+    The deformation equations are based on a 2D elastic sheet with a constant
+    Poisson's ratio. The data can then be predicted at any desired location.
+
+    The east and north data components are coupled through the elastic
+    deformation equations. This coupling is controlled by the Poisson's ratio,
+    which is usually between -1 and 1. The special case of Poisson's ratio -1
+    leads to an uncoupled interpolation, meaning that the east and north
+    components don't interfere with each other.
+
+    The point forces are traditionally placed under each data point. The force
+    locations are set the first time :meth:`~verde.Vector2D.fit` is called.
+    Subsequent calls will fit using the same force locations as the first call.
+    This configuration results in an exact prediction at the data points but
+    can be unstable.
+
+    [SandwellWessel2016]_ stabilize the solution using Singular Value
+    Decomposition but we use ridge regression instead. The regularization can
+    be controlled using the *damping* argument. Alternatively, we also allow
+    forces to be placed on a regular grid using the *spacing*, *shape*, and/or
+    *region* arguments. Regularization or forces on a grid will result in a
+    least-squares estimate at the data points, not an exact solution. Note that
+    the least-squares solution is required for data weights to have any effect.
+
+    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
+    ----------
+    poisson : float
+        The Poisson's ratio for the elastic deformation Green's functions.
+        Default is 0.5. A value of -1 will lead to uncoupled interpolation of
+        the east and north data components.
+    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.
+    region : None or tuple
+        If not None, then the boundaries (``[W, E, S, N]``) used to generate a
+        regular grid of forces. If None is given, then will use the boundaries
+        of data given to the first call to :meth:`~verde.Vector2D.fit`.
+
+    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.Vector2D.grid` and :meth:`~verde.Vector2D.scatter`
+        methods.
+
+    See also
+    --------
+    verde.vector2d_jacobian : Jacobian matrix for the 2D elastic deformation
+
+    """
+
+    def __init__(self, poisson=0.5, fudge=1e-5, damping=None, shape=None,
+                 spacing=None, region=None):
+        self.poisson = poisson
+        super().__init__(fudge=fudge, damping=damping, shape=shape,
+                         spacing=spacing, region=region)
+
+    def fit(self, coordinates, data, weights=None):
+        """
+        Fit the gridder to the given 2-component vector data.
+
+        The data region is captured and used as default for the
+        :meth:`~verde.Vector2D.grid` and :meth:`~verde.Vector2D.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 : tuple of array
+            A tuple ``(east_component, north_component)`` of arrays with the
+            vector data values at each point.
+        weights : None or tuple array
+            If not None, then the weights assigned to each data point. Must be
+            one array per data component. 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, unpack=False)
+        if len(data) != 2:
+            raise ValueError("Need two data components. Only {} given."
+                             .format(len(data)))
+        # Capture the data region to use as a default when gridding.
+        self.region_ = get_region(coordinates[:2])
+        self.force_coords_ = self._get_force_coordinates(coordinates)
+        if any(w is not None for w in weights):
+            weights = np.concatenate([i.ravel() for i in weights])
+        else:
+            weights = None
+        self._check_weighted_exact_solution(weights)
+        data = np.concatenate([i.ravel() for i in data])
+        jacobian = vector2d_jacobian(coordinates[:2], self.force_coords_,
+                                     self.poisson, self.fudge)
+        self.force_ = self._estimate_forces(jacobian, data, weights)
+        return self
+
+    def predict(self, coordinates):
+        """
+        Evaluate the fitted gridder on the given set of points.
+
+        Requires a fitted estimator (see :meth:`~verde.Vector2D.fit`).
+
+        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 : tuple of arrays
+            A tuple ``(east_component, north_component)`` of arrays with the
+            predicted vector data values at each point.
+
+        """
+        check_is_fitted(self, ['force_', 'force_coords_'])
+        jac = vector2d_jacobian(coordinates[:2], self.force_coords_,
+                                self.poisson, self.fudge)
+        cast = np.broadcast(*coordinates[:2])
+        npoints = cast.size
+        components = jac.dot(self.force_).reshape((2, npoints))
+        return tuple(comp.reshape(cast.shape) for comp in components)
+
+
+def vector2d_jacobian(coordinates, force_coordinates, poisson, fudge=1e-5,
+                      dtype="float32"):
+    """
+    Make the Jacobian matrix for the 2D coupled elastic deformation.
+
+    Follows [SandwellWessel2016]_.
+
+    The Jacobian is segmented into 4 parts, each relating a force component to
+    a data component::
+
+        | J_ee  J_ne |*|f_e| = |d_e|
+        | J_ne  J_nn | |f_n|   |d_n|
+
+    The forces and data are assumed to be stacked into 1D arrays with the east
+    component on top of the north component.
+
+    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.
+    force_coordinates : tuple of arrays
+        Arrays with the coordinates of each vertical force. Should be in the
+        following order: (easting, northing, vertical, ...). Only easting and
+        northing will be used, all subsequent coordinates will be ignored.
+    poisson ; float
+        The Poisson's ratio for the elastic deformation Green's functions.
+        A value of -1 will lead to uncoupled interpolation of
+        the east and north data components (the ``J_ne`` component of the
+        Jacobian is null).
+    fudge : float
+        The positive fudge factor applied to the Green's function to avoid
+        singularities.
+    dtype : str or numpy dtype
+        The type of the Jacobian array.
+
+    Returns
+    -------
+    jacobian : 2D array
+        The (n_data*2, n_forces*2) Jacobian matrix.
+
+    See also
+    --------
+    verde.Vector2D : Coupled gridder for 2-component vector data
+
+    """
+    force_coordinates = [np.atleast_1d(i).ravel()
+                         for i in force_coordinates[:2]]
+    coordinates = [np.atleast_1d(i).ravel() for i in coordinates[:2]]
+    npoints = coordinates[0].size
+    nforces = force_coordinates[0].size
+    # Reshaping the data coordinates to a column vector will automatically
+    # build a distance matrix between each data point and force.
+    east, north = (datac.reshape((npoints, 1)) - forcec
+                   for datac, forcec in zip(coordinates, force_coordinates))
+    distance = np.hypot(east, north, dtype=dtype)
+    # The fudge factor helps avoid singular matrices when the force and
+    # computation point are too close
+    distance += fudge
+    # Pre-compute common terms for the Green's functions of each component
+    ln_r = (3 - poisson)*np.log(distance)
+    over_r2 = (1 + poisson)/distance**2
+    jac = np.empty((npoints*2, nforces*2), dtype=dtype)
+    jac[:npoints, :nforces] = ln_r + over_r2*north**2  # J_ee
+    jac[npoints:, nforces:] = ln_r + over_r2*east**2  # J_nn
+    jac[:npoints, nforces:] = -over_r2*east*north  # J_ne
+    jac[npoints:, :nforces] = jac[:npoints, nforces:]  # J is symmetric
+    return jac


Test patch:
diff --git a/verde/tests/test_chain.py b/verde/tests/test_chain.py
index 7590d02..f804de9 100644
--- a/verde/tests/test_chain.py
+++ b/verde/tests/test_chain.py
@@ -2,12 +2,15 @@
 Test the Chain class
 """
 import numpy.testing as npt
+import numpy as np
 
 from ..datasets.synthetic import CheckerBoard
+from ..blockreduce import BlockReduce
 from ..chain import Chain
 from ..scipygridder import ScipyGridder
 from ..spline import Spline
-from ..trend import Trend
+from ..trend import Trend, VectorTrend
+from ..coordinates import grid_coordinates
 
 
 def test_chain():
@@ -52,3 +55,18 @@ def test_chain_double():
 
     npt.assert_allclose(grd.predict(coords), all_data)
     npt.assert_allclose(grd.score(coords, all_data), 1)
+
+
+def test_chain_vector():
+    "Test chains with vector data"
+    coords = grid_coordinates((-10, 0, 0, 10), spacing=0.1)
+    coefs = [-50, -0.5, 2.1]
+    trend = coefs[0] + coefs[1]*coords[0] + coefs[2]*coords[1]
+    data = (trend, trend.copy())
+    chain = Chain([('mean', BlockReduce(np.mean, spacing=0.5)),
+                   ('trend', VectorTrend(degree=1))])
+    chain.fit(coords, data)
+    npt.assert_allclose(chain.score(coords, data), 1)
+    npt.assert_allclose(chain.predict(coords), data)
+    npt.assert_allclose(chain.named_steps['trend'].component_[0].coef_, coefs)
+    npt.assert_allclose(chain.named_steps['trend'].component_[1].coef_, coefs)
diff --git a/verde/tests/test_spline.py b/verde/tests/test_spline.py
index 1587ce3..4396ad5 100644
--- a/verde/tests/test_spline.py
+++ b/verde/tests/test_spline.py
@@ -80,14 +80,25 @@ def test_spline_grid():
                         rtol=5e-2)
 
 
-def test_spline_warns():
+def test_spline_warns_weights():
     "Check that a warning is issued when using weights and the exact solution."
     data = CheckerBoard().scatter(random_state=100)
     weights = np.ones_like(data.scalars)
     grd = Spline()
-    msg = "Weights are ignored for the exact spline solution"
+    msg = "Weights are ignored for the exact 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
+
+
+def test_spline_warns_underdetermined():
+    "Check that a warning is issued when the problem is underdetermined"
+    data = CheckerBoard().scatter(size=50, random_state=100)
+    grd = Spline(shape=(10, 10))
+    with warnings.catch_warnings(record=True) as warn:
+        grd.fit((data.easting, data.northing), data.scalars)
+        assert len(warn) == 1
+        assert issubclass(warn[-1].category, UserWarning)
+        assert str(warn[-1].message).startswith("Under-determined problem")
diff --git a/verde/tests/test_vector.py b/verde/tests/test_vector.py
new file mode 100644
index 0000000..f249969
--- /dev/null
+++ b/verde/tests/test_vector.py
@@ -0,0 +1,53 @@
+# pylint: disable=redefined-outer-name
+"""
+Test the vector data interpolators
+"""
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+from ..datasets.synthetic import CheckerBoard
+from ..coordinates import grid_coordinates
+from ..vector import Vector2D
+
+
+@pytest.fixture
+def data2d():
+    "Make 2 component vector data"
+    synth = CheckerBoard(region=(2, 8, -3, 3), w_east=20, w_north=20)
+    coordinates = grid_coordinates(synth.region, spacing=0.5)
+    data = tuple(synth.predict(coordinates) for i in range(2))
+    return coordinates, data
+
+
+def test_vector2d(data2d):
+    "See if the exact solution works"
+    coords, data = data2d
+    spline = Vector2D().fit(coords, data)
+    npt.assert_allclose(spline.score(coords, data), 1)
+    npt.assert_allclose(spline.predict(coords), data, rtol=1e-3)
+    # There should be 1 force per data point
+    assert data[0].size == spline.force_coords_[0].size
+    assert data[0].size*2 == spline.force_.size
+    npt.assert_allclose(spline.force_coords_, coords)
+
+
+def test_vector2d_weights(data2d):
+    "Use unit weights and a regular grid solution"
+    coords, data = data2d
+    weights = tuple(np.ones_like(data[0]) for i in range(2))
+    spline = Vector2D(shape=(11, 11)).fit(coords, data, weights)
+    npt.assert_allclose(spline.score(coords, data), 1, rtol=1e-3)
+    npt.assert_allclose(spline.predict(coords), data, rtol=5e-2)
+    assert spline.force_coords_[0].size == 11*11
+    assert spline.force_.size == 2*11*11
+
+
+def test_vector2d_fails(data2d):
+    "Should fail if not given 2 data components"
+    coords, data = data2d
+    spline = Vector2D()
+    with pytest.raises(ValueError):
+        spline.fit(coords, data[0])
+    with pytest.raises(ValueError):
+        spline.fit(coords, data + coords)