Repository: fatiando/verde
Pull request #56 (2018-05-22 03:28:32)

Standardize functions to use coordinates tuples
Some functions where using `easting, northing` as input. I changed all
of them to take in a `coordinates` tuple instead. Now everything takes
it and uses the first 2 elements. Return values are also grouped into 
coordinates tuples.


Patch:
diff --git a/examples/block_reduce.py b/examples/block_reduce.py
index c1421fd..316a661 100644
--- a/examples/block_reduce.py
+++ b/examples/block_reduce.py
@@ -19,9 +19,10 @@ import verde as vd
 data = vd.datasets.fetch_baja_bathymetry()
 
 # Decimate the data using a blocked median with 10 arc-minute blocks
-lon, lat, bathymetry = vd.block_reduce(
-    data.longitude, data.latitude, data.bathymetry_m,
+coordinates, bathymetry = vd.block_reduce(
+    coordinates=(data.longitude, data.latitude), data=data.bathymetry_m,
     reduction=np.median, spacing=10/60)
+lon, lat = coordinates
 
 print("Original data size:", data.bathymetry_m.size)
 print("Decimated data size:", bathymetry.size)
diff --git a/examples/chain_trend.py b/examples/chain_trend.py
index 8e42cb3..415436d 100644
--- a/examples/chain_trend.py
+++ b/examples/chain_trend.py
@@ -19,19 +19,19 @@ import verde as vd
 
 # Load the Rio de Janeiro total field magnetic anomaly data
 data = vd.datasets.fetch_rio_magnetic_anomaly()
-region = vd.get_region(data.longitude, data.latitude)
+region = vd.get_region((data.longitude, data.latitude))
 
 # Before gridding, we need to decimate the data using a blocked mean to avoid
 # aliasing
 spacing = 0.5/60
-lon, lat, total_field = vd.block_reduce(data.longitude, data.latitude,
-                                        data.total_field_anomaly_nt,
-                                        reduction=np.median, spacing=spacing)
+coordinates, total_field = vd.block_reduce(
+    (data.longitude, data.latitude), data.total_field_anomaly_nt,
+    reduction=np.median, spacing=spacing)
 
 # 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)
+proj_coordinates = projection(*coordinates)
 
 # Create a chain that fits a 2nd degree trend to the anomaly data and then a
 # standard gridder to the residuals
@@ -39,7 +39,7 @@ chain = vd.Chain([('trend', vd.Trend(degree=2)),
                   ('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)
+chain.fit(proj_coordinates, total_field)
 
 # Each component of the chain can be accessed separately using the
 # 'named_steps' attribute
@@ -73,7 +73,7 @@ pc = ax.pcolormesh(grid.longitude, grid.latitude, grid.total_field_anomaly,
 cb = plt.colorbar(pc, pad=0.01)
 cb.set_label('total field anomaly [nT]')
 # Plot the locations of the decimated data
-ax.plot(lon, lat, '.k', markersize=0.5, transform=crs)
+ax.plot(*coordinates, '.k', markersize=0.5, transform=crs)
 # Set the proper ticks for a Cartopy map
 ax.set_xticks(np.arange(-42.5, -42, 0.1), crs=crs)
 ax.set_yticks(np.arange(-22.4, -22, 0.1), crs=crs)
diff --git a/examples/data/sample_california_gps.py b/examples/data/sample_california_gps.py
index 3714f9b..2c0547f 100644
--- a/examples/data/sample_california_gps.py
+++ b/examples/data/sample_california_gps.py
@@ -28,7 +28,7 @@ def setup_map(ax):
     ax.set_yticks(np.arange(33, 42, 2), crs=crs)
     ax.xaxis.set_major_formatter(LongitudeFormatter())
     ax.yaxis.set_major_formatter(LatitudeFormatter())
-    ax.set_extent(vd.get_region(data.longitude, data.latitude), crs=crs)
+    ax.set_extent(vd.get_region((data.longitude, data.latitude)), crs=crs)
     # Plot the land and ocean as a solid color
     ax.add_feature(cfeature.LAND)
     ax.add_feature(cfeature.OCEAN)
diff --git a/examples/data/sample_rio_magnetic_anomaly.py b/examples/data/sample_rio_magnetic_anomaly.py
index 376e6ca..8453c13 100644
--- a/examples/data/sample_rio_magnetic_anomaly.py
+++ b/examples/data/sample_rio_magnetic_anomaly.py
@@ -41,6 +41,6 @@ ax.set_yticks(np.arange(-22.4, -22, 0.1), crs=crs)
 ax.xaxis.set_major_formatter(LongitudeFormatter())
 ax.yaxis.set_major_formatter(LatitudeFormatter())
 # Set the extent of the plot to the limits of the data
-ax.set_extent(vd.get_region(data.longitude, data.latitude), crs=crs)
+ax.set_extent(vd.get_region((data.longitude, data.latitude)), crs=crs)
 plt.tight_layout()
 plt.show()
diff --git a/examples/distance_mask.py b/examples/distance_mask.py
index 2a701ab..e67ca0c 100644
--- a/examples/distance_mask.py
+++ b/examples/distance_mask.py
@@ -15,17 +15,17 @@ import verde as vd
 # The Baja California bathymetry dataset has big gaps on land. We want to mask
 # these gaps on a grid that we'll generate over the region
 data = vd.datasets.fetch_baja_bathymetry()
-region = vd.get_region(data.longitude, data.latitude)
+region = vd.get_region((data.longitude, data.latitude))
 
 # Generate the coordinates and a dummy grid with only ones that we'll want to
 # mask
 spacing = 10/60
-lons, lats = vd.grid_coordinates(region, spacing=spacing)
-dummy_data = np.ones_like(lons)
+coordinates = vd.grid_coordinates(region, spacing=spacing)
+dummy_data = np.ones_like(coordinates[0])
 
 # Generate a mask for points that are more than 2 grid spacings away from any
 # data point. The mask is True for points that are too far.
-mask = vd.distance_mask(lons, lats, data.longitude, data.latitude,
+mask = vd.distance_mask(coordinates, (data.longitude, data.latitude),
                         maxdist=spacing*2)
 print(mask)
 
@@ -38,7 +38,7 @@ plt.figure(figsize=(7, 6))
 ax = plt.axes(projection=ccrs.Mercator())
 ax.set_title('Only grid points that are close to data', pad=25)
 ax.plot(data.longitude, data.latitude, '.y', markersize=0.5, transform=crs)
-ax.pcolormesh(lons, lats, dummy_data, transform=crs)
+ax.pcolormesh(*coordinates, dummy_data, transform=crs)
 ax.gridlines(draw_labels=True)
 ax.set_extent(region, crs=crs)
 plt.tight_layout()
diff --git a/examples/scipygridder.py b/examples/scipygridder.py
index d04e8b2..a285608 100644
--- a/examples/scipygridder.py
+++ b/examples/scipygridder.py
@@ -31,21 +31,21 @@ data = vd.datasets.fetch_baja_bathymetry()
 # the oversampling along the ship tracks. We'll use a blocked median with 5
 # arc-minute blocks.
 spacing = 5/60
-lon, lat, bathymetry = vd.block_reduce(data.longitude, data.latitude,
-                                       data.bathymetry_m, reduction=np.median,
-                                       spacing=spacing)
+coordinates, bathymetry = vd.block_reduce(
+    (data.longitude, data.latitude), data.bathymetry_m, reduction=np.median,
+    spacing=spacing)
 
 # 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)
+proj_coordinates = projection(*coordinates)
 
 # Now we can set up a gridder for the decimated data
-grd = vd.ScipyGridder(method='cubic').fit(coordinates, bathymetry)
+grd = vd.ScipyGridder(method='cubic').fit(proj_coordinates, bathymetry)
 print("Gridder used:", grd)
 
 # Get the grid region in geographic coordinates
-region = vd.get_region(data.longitude, data.latitude)
+region = vd.get_region((data.longitude, data.latitude))
 print("Data region:", region)
 
 # The 'grid' method can still make a geographic grid if we pass in a projection
@@ -75,7 +75,7 @@ pc = ax.pcolormesh(grid.longitude, grid.latitude, grid.bathymetry_m,
 cb = plt.colorbar(pc, pad=0.08)
 cb.set_label('bathymetry [m]')
 # Plot the locations of the decimated data
-ax.plot(lon, lat, '.k', markersize=0.5, transform=crs)
+ax.plot(*coordinates, '.k', markersize=0.5, transform=crs)
 # Plot the land as a solid color
 ax.add_feature(cfeature.LAND, edgecolor='black', zorder=2)
 ax.set_extent(region, crs=crs)
diff --git a/examples/spline.py b/examples/spline.py
index f08da01..f75d0fc 100644
--- a/examples/spline.py
+++ b/examples/spline.py
@@ -21,13 +21,13 @@ 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)
+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)
+coordinates, 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)
@@ -37,7 +37,7 @@ 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)
+proj_coordinates = projection(*coordinates)
 
 # Create an array of weights and set the weights for the outliers to a very low
 # value
@@ -46,12 +46,13 @@ 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)
+spline = vd.Spline().fit(proj_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)
+spline_weights = vd.Spline(damping=1e-8).fit(proj_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,
@@ -72,6 +73,7 @@ for ax, title, grid in zip(axes, titles, grids):
     # Plot the land as a solid color
     ax.add_feature(cfeature.LAND, edgecolor='black', zorder=2)
     # Plot the locations of the decimated data
+    lon, lat = coordinates
     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,
diff --git a/examples/trend.py b/examples/trend.py
index 70d104a..c703f25 100644
--- a/examples/trend.py
+++ b/examples/trend.py
@@ -51,7 +51,7 @@ def plot_data(column, i, title):
     ax.xaxis.set_major_formatter(LongitudeFormatter())
     ax.yaxis.set_major_formatter(LatitudeFormatter())
     # Set the plot region to be tight around the data
-    ax.set_extent(vd.get_region(data.longitude, data.latitude))
+    ax.set_extent(vd.get_region((data.longitude, data.latitude)))
     return mappable
 
 
diff --git a/examples/vector_trend.py b/examples/vector_trend.py
index c4e0110..f8edd72 100644
--- a/examples/vector_trend.py
+++ b/examples/vector_trend.py
@@ -69,7 +69,7 @@ def plot_trend(ax, component, title):
     ax.xaxis.set_major_formatter(LongitudeFormatter())
     ax.yaxis.set_major_formatter(LatitudeFormatter())
     # ax.coastlines(color='white')
-    ax.set_extent(vd.get_region(data.longitude, data.latitude), crs=crs)
+    ax.set_extent(vd.get_region((data.longitude, data.latitude)), crs=crs)
 
 
 # Make a plot of the data using Cartopy to handle projections and coastlines
diff --git a/verde/base.py b/verde/base.py
index 44b1fb2..03179ef 100644
--- a/verde/base.py
+++ b/verde/base.py
@@ -219,7 +219,7 @@ class BaseGridder(BaseEstimator):
         if shape is None and spacing is None:
             shape = (101, 101)
         dims = get_dims(self, dims)
-        region = get_region(self, region)
+        region = get_instance_region(self, region)
         coordinates = grid_coordinates(region, shape=shape, spacing=spacing,
                                        adjust=adjust)
         if projection is None:
@@ -276,7 +276,7 @@ class BaseGridder(BaseEstimator):
 
         """
         dims = get_dims(self, dims)
-        region = get_region(self, region)
+        region = get_instance_region(self, region)
         coordinates = scatter_points(region, size, random_state)
         if projection is None:
             data = check_data(self.predict(coordinates))
@@ -410,7 +410,7 @@ def get_data_names(data, data_names):
     return data_types[len(data) - 1]
 
 
-def get_region(instance, region):
+def get_instance_region(instance, region):
     """
     Get the region attribute stored in instance if one is not provided.
     """
diff --git a/verde/chain.py b/verde/chain.py
index 4b4b5fa..fb04246 100644
--- a/verde/chain.py
+++ b/verde/chain.py
@@ -79,7 +79,7 @@ class Chain(BaseGridder):
             Returns this estimator instance for chaining operations.
 
         """
-        self.region_ = get_region(*coordinates[:2])
+        self.region_ = get_region(coordinates[:2])
         residuals = data
         for _, step in self.steps:
             step.fit(coordinates, residuals, weights)
diff --git a/verde/coordinates.py b/verde/coordinates.py
index 0dd3c1c..893a8bd 100644
--- a/verde/coordinates.py
+++ b/verde/coordinates.py
@@ -36,16 +36,16 @@ def check_region(region):
                          .format(region))
 
 
-def get_region(easting, northing):
+def get_region(coordinates):
     """
     Get the bounding region of the given coordinates.
 
     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.
+    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
     -------
@@ -56,11 +56,12 @@ def get_region(easting, northing):
     Examples
     --------
 
-    >>> east, north = grid_coordinates((0, 1, -10, -6), shape=(10, 10))
-    >>> print(get_region(east, north))
+    >>> coords = grid_coordinates((0, 1, -10, -6), shape=(10, 10))
+    >>> print(get_region(coords))
     (0.0, 1.0, -10.0, -6.0)
 
     """
+    easting, northing = coordinates[:2]
     region = (np.min(easting), np.max(easting),
               np.min(northing), np.max(northing))
     return region
@@ -423,7 +424,7 @@ def profile_coordinates(point1, point2, size, coordinate_system='cartesian'):
     return easting, northing, distances
 
 
-def inside(easting, northing, region, out=None, tmp=None):
+def inside(coordinates, region, out=None, tmp=None):
     """
     Determine which points fall inside a given region.
 
@@ -431,10 +432,10 @@ def inside(easting, northing, region, out=None, tmp=None):
 
     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.
+    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.
     region : list = [W, E, S, N]
         The boundaries of a given region in Cartesian or geographic
         coordinates.
@@ -462,7 +463,7 @@ def inside(easting, northing, region, out=None, tmp=None):
     >>> east = np.array([1, 2, 3, 4, 5, 6])
     >>> north = np.array([10, 11, 12, 13, 14, 15])
     >>> region = [2.5, 5.5, 12, 15]
-    >>> print(inside(east, north, region))
+    >>> print(inside((east, north), region))
     [False False  True  True  True False]
     >>> # This also works for 2D-arrays
     >>> east = np.array([[1, 1, 1],
@@ -472,7 +473,7 @@ def inside(easting, northing, region, out=None, tmp=None):
     ...                   [5, 7, 9],
     ...                   [5, 7, 9]])
     >>> region = [0.5, 2.5, 6, 9]
-    >>> print(inside(east, north, region))
+    >>> print(inside((east, north), region))
     [[False  True  True]
      [False  True  True]
      [False False False]]
@@ -480,6 +481,7 @@ def inside(easting, northing, region, out=None, tmp=None):
     """
     check_region(region)
     w, e, s, n = region
+    easting, northing = coordinates[:2]
     # Allocate temporary arrays to minimize memory allocation overhead
     if out is None:
         out = np.empty_like(easting, dtype=np.bool)
diff --git a/verde/grid_math.py b/verde/grid_math.py
index aa5e9db..39a469f 100644
--- a/verde/grid_math.py
+++ b/verde/grid_math.py
@@ -68,7 +68,7 @@ def block_split(coordinates, spacing, adjust='spacing', region=None):
     """
     easting, northing = coordinates[:2]
     if region is None:
-        region = get_region(easting, northing)
+        region = get_region((easting, northing))
     block_coords = tuple(
         i.ravel() for i in grid_coordinates(
             region, spacing=spacing, adjust=adjust, pixel_register=True)
@@ -79,7 +79,7 @@ def block_split(coordinates, spacing, adjust='spacing', region=None):
     return block_coords, labels
 
 
-def block_reduce(easting, northing, data, reduction, spacing, region=None,
+def block_reduce(coordinates, data, reduction, spacing, region=None,
                  adjust='spacing', center_coordinates=False):
     """
     Apply a reduction operation to the data in blocks (windows).
@@ -102,10 +102,10 @@ def block_reduce(easting, northing, data, reduction, spacing, region=None,
 
     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.
+    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 at each point.
     reduction : function
@@ -129,33 +129,36 @@ def block_reduce(easting, northing, data, reduction, spacing, region=None,
 
     Returns
     -------
-    easting, northing, data : arrays
-        The reduced coordinates and data values.
+    blocked_coordinates : tuple of arrays
+        (easting, northing) arrays with the coordinates of each block that
+        contains data.
+    blocked_data : array
+        The block reduced data values.
 
     See also
     --------
-    block_inside : Split a region into blocks and label points accordingly.
+    block_split : Split a region into blocks and label points accordingly.
 
     """
+    easting, northing = coordinates[:2]
     block_coords, labels = block_split((easting, northing), spacing, adjust,
                                        region)
     if center_coordinates:
         table = pd.DataFrame(dict(data=data.ravel(), block=labels))
         blocked = table.groupby('block').aggregate(reduction)
         unique = np.unique(labels)
-        block_east, block_north = [i[unique] for i in block_coords]
+        blocked_coords = tuple(i[unique] for i in block_coords)
     else:
         table = pd.DataFrame(dict(easting=easting.ravel(),
                                   northing=northing.ravel(),
                                   data=data.ravel(),
                                   block=labels))
         blocked = table.groupby('block').aggregate(reduction)
-        block_east = blocked.easting.values
-        block_north = blocked.northing.values
-    return block_east, block_north, blocked.data.values
+        blocked_coords = (blocked.easting.values, blocked.northing.values)
+    return blocked_coords, blocked.data.values
 
 
-def distance_mask(easting, northing, data_easting, data_northing, maxdist):
+def distance_mask(coordinates, data_coordinates, maxdist):
     """
     Create a mask for points that are too far from the given data points.
 
@@ -164,14 +167,13 @@ def distance_mask(easting, northing, data_easting, data_northing, maxdist):
 
     Parameters
     ----------
-    easting : array
-        The West-East coordinates of the points that will be masked.
-    northing : array
-        The South-North coordinates of the points that will be masked.
-    data_easting : array
-        The West-East coordinates of the data points.
-    data_northing : array
-        The South-North coordinates of the data points.
+    coordinates : tuple of arrays
+        Arrays with the coordinates of each point that will be masked. Should
+        be in the following order: (easting, northing, vertical, ...). Only
+        easting and northing will be used, all subsequent coordinates will be
+        ignored.
+    data_coordinates : tuple of arrays
+        Same as *coordinates* but for the data points.
     maxdist : float
         The maximum distance that a point can be from the closest data point.
 
@@ -184,8 +186,8 @@ def distance_mask(easting, northing, data_easting, data_northing, maxdist):
     --------
 
     >>> from verde import grid_coordinates
-    >>> east, north = grid_coordinates((0, 5, -10, -5), spacing=1)
-    >>> mask = distance_mask(east, north, 2.5, -7.5, maxdist=2)
+    >>> coords = grid_coordinates((0, 5, -10, -5), spacing=1)
+    >>> mask = distance_mask(coords, (2.5, -7.5), maxdist=2)
     >>> print(mask)
     [[ True  True  True  True  True  True]
      [ True  True False False  True  True]
@@ -195,6 +197,8 @@ def distance_mask(easting, northing, data_easting, data_northing, maxdist):
      [ True  True  True  True  True  True]]
 
     """
+    data_easting, data_northing = data_coordinates[:2]
+    easting, northing = coordinates[:2]
     data_easting = np.atleast_1d(data_easting)
     data_northing = np.atleast_1d(data_northing)
     data_points = np.transpose((data_easting.ravel(), data_northing.ravel()))
diff --git a/verde/scipy_bridge.py b/verde/scipy_bridge.py
index d8e4840..1f64359 100644
--- a/verde/scipy_bridge.py
+++ b/verde/scipy_bridge.py
@@ -98,7 +98,7 @@ class ScipyGridder(BaseGridder):
         coordinates, data, weights = check_fit_input(coordinates, data,
                                                      weights)
         easting, northing = coordinates[:2]
-        self.region_ = get_region(easting, northing)
+        self.region_ = get_region((easting, northing))
         points = np.column_stack((np.ravel(easting), np.ravel(northing)))
         self.interpolator_ = classes[self.method](points, np.ravel(data),
                                                   **kwargs)
diff --git a/verde/spline.py b/verde/spline.py
index 270e0c0..4f9850b 100644
--- a/verde/spline.py
+++ b/verde/spline.py
@@ -128,17 +128,16 @@ class Spline(BaseGridder):
             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)
+        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_ = (easting.copy(), northing.copy())
+            self.force_coords_ = tuple(i.copy() for i in coordinates[:2])
         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)
+            self.force_coords_ = grid_coordinates(self.region_,
+                                                  shape=self.shape,
+                                                  spacing=self.spacing)
+        jac = spline_jacobian(coordinates[:2], self.force_coords_, self.fudge)
         scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
         jac = scaler.fit_transform(jac)
         if self.damping is None:
@@ -171,17 +170,13 @@ class Spline(BaseGridder):
             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
+        jac = spline_jacobian(coordinates[:2], self.force_coords_, self.fudge)
+        shape = np.broadcast(*coordinates[:2]).shape
         return jac.dot(self.force_).reshape(shape)
 
 
-def spline_jacobian(easting, northing, force_easting, force_northing,
-                    fudge=1e-5):
+def spline_jacobian(coordinates, force_coordinates, fudge=1e-5):
     """
     Make the Jacobian matrix for the 2D biharmonic spline.
 
@@ -190,14 +185,14 @@ def spline_jacobian(easting, northing, force_easting, force_northing,
 
     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.
+    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.
     fudge : float
         The positive fudge factor applied to the Green's function to avoid
         singularities.
@@ -212,8 +207,10 @@ def spline_jacobian(easting, northing, force_easting, force_northing,
     verde.Spline : Biharmonic spline gridder
 
     """
-    easting = np.atleast_1d(easting)
-    northing = np.atleast_1d(northing)
+    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)
     # Reshaping the data to a column vector will automatically build a
diff --git a/verde/trend.py b/verde/trend.py
index f9f4897..bf5c6b1 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -103,9 +103,8 @@ class Trend(BaseGridder):
         """
         coordinates, data, weights = check_fit_input(coordinates, data,
                                                      weights)
-        easting, northing = coordinates[:2]
-        self.region_ = get_region(easting, northing)
-        jac = trend_jacobian(easting, northing, degree=self.degree,
+        self.region_ = get_region(coordinates[:2])
+        jac = trend_jacobian(coordinates[:2], degree=self.degree,
                              dtype=data.dtype)
         scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
         jac = scaler.fit_transform(jac)
@@ -136,9 +135,8 @@ class Trend(BaseGridder):
 
         """
         check_is_fitted(self, ['coef_'])
-        easting, northing = coordinates[:2]
-        jac = trend_jacobian(easting, northing, degree=self.degree)
-        shape = np.broadcast(easting, northing).shape
+        jac = trend_jacobian(coordinates[:2], degree=self.degree)
+        shape = np.broadcast(*coordinates[:2]).shape
         return jac.dot(self.coef_).reshape(shape)
 
 
@@ -179,7 +177,7 @@ def polynomial_power_combinations(degree):
     return tuple(sorted(combinations, key=sum))
 
 
-def trend_jacobian(easting, northing, degree, dtype='float64'):
+def trend_jacobian(coordinates, degree, dtype='float64'):
     """
     Make the Jacobian for a 2D polynomial of the given degree.
 
@@ -188,10 +186,10 @@ def trend_jacobian(easting, northing, degree, dtype='float64'):
 
     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.
+    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.
     degree : int
         The degree of the 2D polynomial. Must be >= 1.
     dtype : str or numpy dtype
@@ -208,13 +206,13 @@ def trend_jacobian(easting, northing, degree, dtype='float64'):
     >>> import numpy as np
     >>> east = np.linspace(0, 4, 5)
     >>> north = np.linspace(-5, -1, 5)
-    >>> print(trend_jacobian(east, north, degree=1, dtype=np.int))
+    >>> print(trend_jacobian((east, north), degree=1, dtype=np.int))
     [[ 1  0 -5]
      [ 1  1 -4]
      [ 1  2 -3]
      [ 1  3 -2]
      [ 1  4 -1]]
-    >>> print(trend_jacobian(east, north, degree=2, dtype=np.int))
+    >>> print(trend_jacobian((east, north), degree=2, dtype=np.int))
     [[ 1  0 -5  0  0 25]
      [ 1  1 -4  1 -4 16]
      [ 1  2 -3  4 -6  9]
@@ -227,6 +225,7 @@ def trend_jacobian(easting, northing, degree, dtype='float64'):
     verde.VectorTrend : Polynomial trend estimator for multi-component data
 
     """
+    easting, northing = coordinates[:2]
     if easting.shape != northing.shape:
         raise ValueError("Coordinate arrays must have the same shape.")
     combinations = polynomial_power_combinations(degree)
@@ -314,7 +313,7 @@ class VectorTrend(BaseGridder):
                 .format(type(weights)))
         coordinates, data, weights = check_fit_input(coordinates, data,
                                                      weights)
-        self.region_ = get_region(*coordinates[:2])
+        self.region_ = get_region(coordinates[:2])
         self.component_ = [
             Trend(degree=self.degree).fit(coordinates, data_comp, weight_comp)
             for data_comp, weight_comp in zip(data, weights)]


Test patch:
diff --git a/verde/tests/test_base.py b/verde/tests/test_base.py
index a21bd18..64a1ddc 100644
--- a/verde/tests/test_base.py
+++ b/verde/tests/test_base.py
@@ -6,8 +6,8 @@ import pytest
 import numpy as np
 import numpy.testing as npt
 
-from ..base import get_dims, get_data_names, get_region, BaseGridder, \
-    check_fit_input
+from ..base import get_dims, get_data_names, get_instance_region, \
+    BaseGridder, check_fit_input
 from ..coordinates import grid_coordinates, scatter_points
 
 
@@ -70,15 +70,15 @@ def test_get_data_names_fails():
         get_data_names(tuple([np.arange(5)]*2), data_names=("meh",))
 
 
-def test_get_region():
-    "Check if get_region finds the correct region"
+def test_get_instance_region():
+    "Check if get_instance_region finds the correct region"
     grd = BaseGridder()
-    assert get_region(grd, region=(1, 2, 3, 4)) == (1, 2, 3, 4)
+    assert get_instance_region(grd, region=(1, 2, 3, 4)) == (1, 2, 3, 4)
     with pytest.raises(ValueError):
-        get_region(grd, region=None)
+        get_instance_region(grd, region=None)
     grd.region_ = (5, 6, 7, 8)
-    assert get_region(grd, region=None) == (5, 6, 7, 8)
-    assert get_region(grd, region=(1, 2, 3, 4)) == (1, 2, 3, 4)
+    assert get_instance_region(grd, region=None) == (5, 6, 7, 8)
+    assert get_instance_region(grd, region=(1, 2, 3, 4)) == (1, 2, 3, 4)
 
 
 class PolyGridder(BaseGridder):
diff --git a/verde/tests/test_grid_math.py b/verde/tests/test_grid_math.py
index 26a25ca..1d43066 100644
--- a/verde/tests/test_grid_math.py
+++ b/verde/tests/test_grid_math.py
@@ -13,12 +13,14 @@ def test_block_reduce():
     region = (-5, 0, 5, 10)
     east, north = grid_coordinates(region, spacing=0.1, pixel_register=True)
     data = 20*np.ones_like(east)
-    re_east, re_north, re_data = block_reduce(
-        east, north, data, np.mean, spacing=1)
-    assert len(re_east) == len(re_north) == len(re_data) == 25
-    npt.assert_allclose(re_data, 20)
-    npt.assert_allclose(re_east[:5], np.linspace(-4.5, -0.5, 5))
-    npt.assert_allclose(re_north[::5], np.linspace(5.5, 9.5, 5))
+    block_coords, block_data = block_reduce((east, north), data, np.mean,
+                                            spacing=1)
+    assert len(block_coords[0]) == 25
+    assert len(block_coords[1]) == 25
+    assert len(block_data) == 25
+    npt.assert_allclose(block_data, 20)
+    npt.assert_allclose(block_coords[0][:5], np.linspace(-4.5, -0.5, 5))
+    npt.assert_allclose(block_coords[1][::5], np.linspace(5.5, 9.5, 5))
 
 
 def test_block_reduce_scatter():
@@ -26,10 +28,12 @@ def test_block_reduce_scatter():
     region = (-5, 0, 5, 10)
     east, north = scatter_points(region, size=10000, random_state=0)
     data = 20*np.ones_like(east)
-    re_east, re_north, re_data = block_reduce(
-        east, north, data, np.mean, spacing=1, region=region,
-        center_coordinates=True)
-    assert len(re_east) == len(re_north) == len(re_data) == 25
-    npt.assert_allclose(re_data, 20)
-    npt.assert_allclose(re_east[:5], np.linspace(-4.5, -0.5, 5))
-    npt.assert_allclose(re_north[::5], np.linspace(5.5, 9.5, 5))
+    block_coords, block_data = block_reduce((east, north), data, np.mean,
+                                            spacing=1, region=region,
+                                            center_coordinates=True)
+    assert len(block_coords[0]) == 25
+    assert len(block_coords[1]) == 25
+    assert len(block_data) == 25
+    npt.assert_allclose(block_data, 20)
+    npt.assert_allclose(block_coords[0][:5], np.linspace(-4.5, -0.5, 5))
+    npt.assert_allclose(block_coords[1][::5], np.linspace(5.5, 9.5, 5))
diff --git a/verde/tests/test_trend.py b/verde/tests/test_trend.py
index 516a1ca..5374c16 100644
--- a/verde/tests/test_trend.py
+++ b/verde/tests/test_trend.py
@@ -73,7 +73,7 @@ def test_trend_jacobian_fails():
     "Test failing conditions for the trend jacobian builder"
     east, north = np.arange(50), np.arange(30)
     with pytest.raises(ValueError):
-        trend_jacobian(east, north, degree=1)
+        trend_jacobian((east, north), degree=1)
 
 
 def test_vector_trend(simple_2d_model):