Add a convex hull masking function

doc/api/index.rst

@@ -61,6 +61,15 @@ Coordinate Manipulation
     inside
     block_split
 
+Masking
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+    distance_mask
+    convexhull_mask
+
 Utilities
 ---------
 
@@ -69,7 +78,6 @@ Utilities
 
     test
     maxabs
-    distance_mask
     variance_to_weights
     grid_to_table
     median_distance

examples/convex_hull_mask.py

@@ -0,0 +1,45 @@
+"""
+Mask grid points by convex hull
+===============================
+
+Sometimes, data points are unevenly distributed. In such cases, we might not
+want to have interpolated grid points that are too far from any data point.
+Function :func:`verde.convexhull_mask` allows us to set grid points that fall
+outside of the convex hull of the data points to NaN or some other value.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import pyproj
+import numpy as np
+import verde as vd
+
+# The Baja California bathymetry dataset has big gaps on land. We want to mask
+# these gaps on a dummy grid that we'll generate over the region.
+data = vd.datasets.fetch_baja_bathymetry()
+region = vd.get_region((data.longitude, data.latitude))
+
+# Generate the coordinates for a regular grid mask
+spacing = 10 / 60
+coordinates = vd.grid_coordinates(region, spacing=spacing)
+
+# Generate a mask for points. The mask is True for points that are within the
+# convex hull.
+mask = vd.convexhull_mask((data.longitude, data.latitude), coordinates=coordinates)
+print(mask)
+
+# Create a dummy grid with ones that we can mask to show the results. Turn
+# points that are outside of the convex hull into NaNs so they won't show up in
+# our plot.
+dummy_data = np.ones_like(coordinates[0])
+dummy_data[~mask] = np.nan
+
+# Make a plot of the masked data and the data locations.
+crs = ccrs.PlateCarree()
+plt.figure(figsize=(7, 6))
+ax = plt.axes(projection=ccrs.Mercator())
+ax.set_title("Only keep grid points that inside of the convex hull")
+ax.plot(data.longitude, data.latitude, ".y", markersize=0.5, transform=crs)
+ax.pcolormesh(*coordinates, dummy_data, transform=crs)
+vd.datasets.setup_baja_bathymetry_map(ax, land=None)
+plt.tight_layout()
+plt.show()

examples/spline_weights.py

@@ -62,12 +62,7 @@
     dims=["latitude", "longitude"],
     data_names=["velocity"],
 )
-grid = vd.distance_mask(
-    (data.longitude, data.latitude),
-    maxdist=5 * spacing * 111e3,
-    grid=grid_full,
-    projection=projection,
-)
+grid = vd.convexhull_mask((data.longitude, data.latitude), grid=grid_full)
 
 fig, axes = plt.subplots(
     1, 2, figsize=(9, 7), subplot_kw=dict(projection=ccrs.Mercator())

examples/vector_uncoupled.py

@@ -61,13 +61,11 @@
 print("Cross-validation R^2 score: {:.2f}".format(score))
 
 # Interpolate the wind speed onto a regular geographic grid and mask the data that are
-# far from the observation points
+# outside of the convex hull of the data points.
 grid_full = chain.grid(
     region, spacing=spacing, projection=projection, dims=["latitude", "longitude"]
 )
-grid = vd.distance_mask(
-    coordinates, maxdist=3 * spacing * 111e3, grid=grid_full, projection=projection
-)
+grid = vd.convexhull_mask(coordinates, grid=grid_full)
 
 # Make maps of the original and gridded wind speed
 plt.figure(figsize=(6, 6))

tutorials/weights.py

@@ -246,6 +246,8 @@ def plot_data(coordinates, velocity, weights, title_data, title_weights):
     dims=["latitude", "longitude"],
     data_names=["velocity"],
 )
+# Avoid showing interpolation too far away data points.
+grid = vd.convexhull_mask(coordinates, grid=grid)
 
 ########################################################################################
 # Calculate an unweighted spline as well for comparison.
@@ -262,6 +264,7 @@ def plot_data(coordinates, velocity, weights, title_data, title_weights):
     dims=["latitude", "longitude"],
     data_names=["velocity"],
 )
+grid_unweighted = vd.convexhull_mask(coordinates, grid=grid_unweighted)
 
 ########################################################################################
 # Finally, plot the weighted and unweighted grids side by side.

verde/__init__.py

@@ -13,7 +13,7 @@
     project_region,
     longitude_continuity,
 )
-from .mask import distance_mask
+from .mask import distance_mask, convexhull_mask
 from .utils import variance_to_weights, maxabs, grid_to_table
 from .io import load_surfer
 from .distances import median_distance

verde/mask.py

@@ -3,6 +3,9 @@
 """
 import numpy as np
 
+# pylint doesn't pick up on this import for some reason
+from scipy.spatial import Delaunay  # pylint: disable=no-name-in-module
+
 from .base import n_1d_arrays
 from .utils import kdtree
 
@@ -93,14 +96,7 @@ def distance_mask(
      [nan nan nan nan nan nan]]
 
     """
-    if coordinates is None and grid is None:
-        raise ValueError("Either coordinates or grid must be given.")
-    if coordinates is None:
-        dims = [grid[var].dims for var in grid.data_vars][0]
-        coordinates = np.meshgrid(grid.coords[dims[1]], grid.coords[dims[0]])
-    if len(set(i.shape for i in coordinates)) != 1:
-        raise ValueError("Coordinate arrays must have the same shape.")
-    shape = coordinates[0].shape
+    coordinates, shape = _get_grid_coordinates(coordinates, grid)
     if projection is not None:
         data_coordinates = projection(*n_1d_arrays(data_coordinates, 2))
         coordinates = projection(*n_1d_arrays(coordinates, 2))
@@ -110,3 +106,100 @@ def distance_mask(
     if grid is not None:
         return grid.where(mask)
     return mask
+
+
+def convexhull_mask(
+    data_coordinates, coordinates=None, grid=None,
+):
+    """
+    Mask grid points that are outside the convex hull of the given data points.
+
+    Either *coordinates* or *grid* must be given:
+
+    * If *coordinates* is not None, produces an array that is False when a
+      point is outside the convex hull and True otherwise.
+    * If *grid* is not None, produces a mask and applies it to *grid* (an
+      :class:`xarray.Dataset`).
+
+    Parameters
+    ----------
+    data_coordinates : tuple of arrays
+        Same as *coordinates* but for the data points.
+    coordinates : None or tuple of arrays
+        Arrays with the coordinates of each point that will be masked. Should
+        be in the following order: (easting, northing, ...). Only easting and
+        northing will be used, all subsequent coordinates will be ignored.
+    grid : None or :class:`xarray.Dataset`
+        2D grid with values to be masked. Will use the first two dimensions of
+        the grid as northing and easting coordinates, respectively. The mask
+        will be applied to *grid* using the :meth:`xarray.Dataset.where`
+        method.
+
+    Returns
+    -------
+    mask : array or :class:`xarray.Dataset`
+        If *coordinates* was given, then a boolean array with the same shape as
+        the elements of *coordinates*. If *grid* was given, then an
+        :class:`xarray.Dataset` with the mask applied to it.
+
+    Examples
+    --------
+
+    >>> from verde import grid_coordinates
+    >>> region = (0, 5, -10, -4)
+    >>> spacing = 1
+    >>> coords = grid_coordinates(region, spacing=spacing)
+    >>> data_coords = ((2, 3, 2, 3), (-9, -9, -6, -6))
+    >>> mask = convexhull_mask(data_coords, coordinates=coords)
+    >>> print(mask)
+    [[False False False False False False]
+     [False False  True  True False False]
+     [False False  True  True False False]
+     [False False  True  True False False]
+     [False False  True  True False False]
+     [False False False False False False]
+     [False False False False False False]]
+    >>> # Mask an xarray.Dataset directly
+    >>> import xarray as xr
+    >>> coords_dict = {"easting": coords[0][0, :], "northing": coords[1][:, 0]}
+    >>> data_vars = {"scalars": (["northing", "easting"], np.ones(mask.shape))}
+    >>> grid = xr.Dataset(data_vars, coords=coords_dict)
+    >>> masked = convexhull_mask(data_coords, grid=grid)
+    >>> print(masked.scalars.values)
+    [[nan nan nan nan nan nan]
+     [nan nan  1.  1. nan nan]
+     [nan nan  1.  1. nan nan]
+     [nan nan  1.  1. nan nan]
+     [nan nan  1.  1. nan nan]
+     [nan nan nan nan nan nan]
+     [nan nan nan nan nan nan]]
+
+    """
+    coordinates, shape = _get_grid_coordinates(coordinates, grid)
+    n_coordinates = len(data_coordinates)
+    triangles = Delaunay(np.transpose(n_1d_arrays(data_coordinates, n_coordinates)))
+    # Find the triangle that contains each grid point.
+    # -1 indicates that it's not in any triangle.
+    in_triangle = triangles.find_simplex(
+        np.transpose(n_1d_arrays(coordinates, n_coordinates))
+    )
+    mask = (in_triangle != -1).reshape(shape)
+    if grid is not None:
+        return grid.where(mask)
+    return mask
+
+
+def _get_grid_coordinates(coordinates, grid):
+    """
+    If coordinates is given, return it and their shape. Otherwise, get
+    coordinate arrays from the grid.
+    """
+    if coordinates is None and grid is None:
+        raise ValueError("Either coordinates or grid must be given.")
+    if coordinates is None:
+        dims = [grid[var].dims for var in grid.data_vars][0]
+        coordinates = np.meshgrid(grid.coords[dims[1]], grid.coords[dims[0]])
+    if len(set(i.shape for i in coordinates)) != 1:
+        raise ValueError("Coordinate arrays must have the same shape.")
+    shape = coordinates[0].shape
+    return coordinates, shape