Add function to project 2D gridded data

doc/api/index.rst

@@ -57,12 +57,20 @@ Coordinate Manipulation
     profile_coordinates
     get_region
     pad_region
-    project_region
     inside
     block_split
     rolling_window
     expanding_window
 
+Projection
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+    project_region
+    project_grid
+
 Masking
 -------
 

examples/project_grid.py

@@ -0,0 +1,78 @@
+"""
+Projection of gridded data
+==========================
+
+Sometimes gridded data products need to be projected before they can be used.
+For example, you might want to project the topography of Antarctica from
+geographic latitude and longitude to a planar polar stereographic projection
+before doing your analysis. When projecting, the data points will likely not
+fall on a regular grid anymore and must be interpolated (re-sampled) onto a new
+grid.
+
+The :func:`verde.project_grid` function automates this process using the
+interpolation methods available in Verde. An input grid
+(:class:`xarray.DataArray`) is interpolated onto a new grid in the given
+`pyproj <https://jswhit.github.io/pyproj/>`__ projection. The function takes
+care of choosing a default grid spacing and region, running a blocked mean to
+avoid spatial aliasing (using :class:`~verde.BlockReduce`), and masking the
+points in the new grid that aren't constrained by the original data (using
+:func:`~verde.convexhull_mask`).
+
+In this example, we'll generate a synthetic geographic grid with a checkerboard
+pattern around the South pole. We'll project the grid to South Polar
+Stereographic, revealing the checkered pattern of the data.
+
+.. note::
+
+    The interpolation methods are limited to what is available in Verde and
+    there is only support for single 2D grids. For more sophisticated use
+    cases, you might want to try
+    `pyresample <https://github.com/pytroll/pyresample>`__ instead.
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import pyproj
+import verde as vd
+
+
+# We'll use synthetic data near the South pole to highlight the effects of the
+# projection. EPSG 3031 is a South Polar Stereographic projection.
+projection = pyproj.Proj("epsg:3031")
+
+# Create a synthetic geographic grid using a checkerboard pattern
+region = (0, 360, -90, -60)
+spacing = 0.25
+wavelength = 10 * 1e5  # The size of the cells in the checkerboard
+checkerboard = vd.datasets.CheckerBoard(
+    region=vd.project_region(region, projection), w_east=wavelength, w_north=wavelength
+)
+data = checkerboard.grid(
+    region=region,
+    spacing=spacing,
+    projection=projection,
+    data_names=["checkerboard"],
+    dims=("latitude", "longitude"),
+)
+print("Geographic grid:")
+print(data)
+
+# Do the projection while setting the output grid spacing (in projected meters). Set
+# the coordinates names to x and y since they aren't really "northing" or
+# "easting".
+polar_data = vd.project_grid(
+    data.checkerboard, projection, spacing=0.5 * 1e5, dims=("y", "x")
+)
+print("\nProjected grid:")
+print(polar_data)
+
+# Plot the original and projected grids
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 6))
+data.checkerboard.plot(
+    ax=ax1, cbar_kwargs=dict(orientation="horizontal", aspect=50, pad=0.1)
+)
+ax1.set_title("Geographic Grid")
+polar_data.plot(ax=ax2, cbar_kwargs=dict(orientation="horizontal", aspect=50, pad=0.1))
+ax2.set_title("Polar Stereographic Grid")
+plt.tight_layout()
+plt.show()

verde/__init__.py

@@ -12,7 +12,6 @@
     profile_coordinates,
     get_region,
     pad_region,
-    project_region,
     longitude_continuity,
 )
 from .mask import distance_mask, convexhull_mask
@@ -26,6 +25,7 @@
 from .spline import Spline, SplineCV
 from .model_selection import cross_val_score, train_test_split
 from .vector import Vector, VectorSpline2D
+from .projections import project_region, project_grid
 
 
 def test(doctest=True, verbose=True, coverage=False, figures=True):

verde/coordinates.py

@@ -110,40 +110,6 @@ def pad_region(region, pad):
     return padded
 
 
-def project_region(region, projection):
-    """
-    Calculate the bounding box of a region in projected coordinates.
-
-    Parameters
-    ----------
-    region : list = [W, E, S, N]
-        The boundaries of a given region in Cartesian or geographic
-        coordinates.
-    projection : callable or None
-        If not None, then should be a callable object (like a function)
-        ``projection(easting, northing) -> (proj_easting, proj_northing)`` that
-        takes in easting and northing coordinate arrays and returns projected
-        northing and easting coordinate arrays.
-
-    Returns
-    -------
-    proj_region : list = [W, E, S, N]
-        The bounding box of the projected region.
-
-    Examples
-    --------
-
-    >>> def projection(x, y):
-    ...     return (2*x, -1*y)
-    >>> project_region((3, 5, -9, -4), projection)
-    (6.0, 10.0, 4.0, 9.0)
-
-    """
-    east, north = grid_coordinates(region, shape=(101, 101))
-    east, north = projection(east.ravel(), north.ravel())
-    return (east.min(), east.max(), north.min(), north.max())
-
-
 def scatter_points(region, size, random_state=None, extra_coords=None):
     """
     Generate the coordinates for a random scatter of points.
@@ -536,6 +502,62 @@ def spacing_to_shape(region, spacing, adjust):
     return (nnorth, neast), (w, e, s, n)
 
 
+def shape_to_spacing(region, shape, pixel_register=False):
+    """
+    Calculate the spacing of a grid given region and shape.
+
+    Parameters
+    ----------
+    region : list = [W, E, S, N]
+        The boundaries of a given region in Cartesian or geographic
+        coordinates.
+    shape : tuple = (n_north, n_east) or None
+        The number of points in the South-North and West-East directions,
+        respectively.
+    pixel_register : bool
+        If True, the coordinates will refer to the center of each grid pixel
+        instead of the grid lines. In practice, this means that there will be
+        one less element per dimension of the grid when compared to grid line
+        registered (only if given *spacing* and not *shape*). Default is False.
+
+    Returns
+    -------
+    spacing : tuple = (s_north, s_east)
+        The grid spacing in the South-North and West-East directions,
+        respectively.
+
+    Examples
+    --------
+
+    >>> spacing = shape_to_spacing([0, 10, -5, 1], (7, 11))
+    >>> print("{:.1f}, {:.1f}".format(*spacing))
+    1.0, 1.0
+    >>> spacing = shape_to_spacing([0, 10, -5, 1], (14, 11))
+    >>> print("{:.1f}, {:.1f}".format(*spacing))
+    0.5, 1.0
+    >>> spacing = shape_to_spacing([0, 10, -5, 1], (7, 21))
+    >>> print("{:.1f}, {:.1f}".format(*spacing))
+    1.0, 0.5
+    >>> spacing = shape_to_spacing(
+    ...     [-0.5, 10.5, -5.5, 1.5], (7, 11), pixel_register=True,
+    ... )
+    >>> print("{:.1f}, {:.1f}".format(*spacing))
+    1.0, 1.0
+    >>> spacing = shape_to_spacing(
+    ...     [-0.25, 10.25, -5.5, 1.5], (7, 21), pixel_register=True,
+    ... )
+    >>> print("{:.1f}, {:.1f}".format(*spacing))
+    1.0, 0.5
+
+    """
+    spacing = []
+    for i, n_points in enumerate(reversed(shape)):
+        if not pixel_register:
+            n_points -= 1
+        spacing.append((region[2 * i + 1] - region[2 * i]) / n_points)
+    return tuple(reversed(spacing))
+
+
 def profile_coordinates(point1, point2, size, extra_coords=None):
     """
     Coordinates for a profile along a straight line between two points.

verde/projections.py

@@ -0,0 +1,161 @@
+"""
+Operations with projections for grids, regions, etc.
+"""
+import numpy as np
+
+from .coordinates import grid_coordinates, get_region, shape_to_spacing, check_region
+from .utils import grid_to_table
+from .scipygridder import ScipyGridder
+from .blockreduce import BlockReduce
+from .chain import Chain
+from .mask import convexhull_mask
+
+
+def project_region(region, projection):
+    """
+    Calculate the bounding box of a region in projected coordinates.
+
+    Parameters
+    ----------
+    region : list = [W, E, S, N]
+        The boundaries of a given region in Cartesian or geographic
+        coordinates.
+    projection : callable
+        Should be a callable object (like a function) ``projection(easting,
+        northing) -> (proj_easting, proj_northing)`` that takes in easting and
+        northing coordinate arrays and returns projected northing and easting
+        coordinate arrays.
+
+    Returns
+    -------
+    proj_region : list = [W, E, S, N]
+        The bounding box of the projected region.
+
+    Examples
+    --------
+
+    >>> def projection(x, y):
+    ...     return (2*x, -1*y)
+    >>> project_region((3, 5, -9, -4), projection)
+    (6.0, 10.0, 4.0, 9.0)
+
+    """
+    east, north = grid_coordinates(region, shape=(101, 101))
+    east, north = projection(east.ravel(), north.ravel())
+    return (east.min(), east.max(), north.min(), north.max())
+
+
+def project_grid(grid, projection, method="linear", antialias=True, **kwargs):
+    """
+    Apply the given map projection to a grid and re-sample it.
+
+    Creates a new grid in the projected coordinates by interpolating the
+    original values using the chosen *method* (linear by default). Before
+    interpolation, apply a blocked mean operation (:class:`~verde.BlockReduce`)
+    to avoid aliasing when the projected coordinates become oversampled in some
+    regions (which would cause the interpolation to down-sample the original
+    data). For example, applying a polar projection results in oversampled data
+    close to the pole.
+
+    Points that fall outside the convex hull of the original data will be
+    masked (see :func:`~verde.convexhull_mask`) since they are not constrained
+    by any data points.
+
+    Any arguments that can be passed to the
+    :meth:`~verde.base.BaseGridder.grid` method of Verde gridders can be passed
+    to this function as well. Use this to set a region and spacing (or shape)
+    for the projected grid. The region and spacing must be in *projected
+    coordinates*.
+
+    If no region is provided, the bounding box of the projected data will be
+    used. If no spacing or shape is provided, the shape of the input *grid*
+    will be used for the projected grid.
+
+    By default, the ``data_names`` argument will be set to the name of the data
+    variable of the input *grid* (if it has been set).
+
+    .. note::
+
+        The interpolation methods are limited to what is available in Verde and
+        there is only support for single 2D grids. For more sophisticated use
+        cases, you might want to try
+        `pyresample <https://github.com/pytroll/pyresample>`__ instead.
+
+    Parameters
+    ----------
+    grid : :class:`xarray.DataArray`
+        A single 2D grid of values. The first dimension is assumed to be the
+        northing/latitude dimension while the second is assumed to be the
+        easting/longitude dimension.
+    projection : callable
+        Should be a callable object (like a function) ``projection(easting,
+        northing) -> (proj_easting, proj_northing)`` that takes in easting and
+        northing coordinate arrays and returns projected northing and easting
+        coordinate arrays.
+    method : string or Verde gridder
+        If a string, will use it to create a :class:`~verde.ScipyGridder` with
+        the corresponding method (nearest, linear, or cubic). Otherwise, should
+        be a gridder/estimator object, like :class:`~verde.Spline`. Default is
+        ``"linear"``.
+    antialias : bool
+        If True, will run a :class:`~verde.BlockReduce` with a mean function to
+        avoid aliasing when the projection results in oversampling of the data
+        in some areas (for example, in polar projections). If False, will not
+        run the blocked mean.
+
+    Returns
+    -------
+    projected_grid : :class:`xarray.DataArray`
+        The projected grid, interpolated with the given parameters.
+
+    """
+    if hasattr(grid, "data_vars"):
+        raise ValueError(
+            "Projecting xarray.Dataset is not currently supported. "
+            "Please provide a DataArray instead."
+        )
+    if len(grid.dims) != 2:
+        raise ValueError(
+            "Projecting grids with number of dimensions other than 2 is not "
+            "currently supported (dimensions of the given DataArray: {}).".format(
+                len(grid.dims)
+            )
+        )
+
+    # Can be set to None for some data arrays depending on how they are created
+    # so we can't just rely on the default value for getattr.
+    name = getattr(grid, "name", None)
+    if name is None:
+        name = "scalars"
+
+    data = grid_to_table(grid).dropna()
+    coordinates = projection(data[grid.dims[1]].values, data[grid.dims[0]].values)
+    data_region = get_region(coordinates)
+
+    region = kwargs.pop("region", data_region)
+    shape = kwargs.pop("shape", grid.shape)
+    spacing = kwargs.pop("spacing", shape_to_spacing(region, shape))
+
+    check_region(region)
+
+    steps = []
+    if antialias:
+        steps.append(
+            ("mean", BlockReduce(np.mean, spacing=spacing, region=data_region))
+        )
+    if isinstance(method, str):
+        steps.append(("spline", ScipyGridder(method)))
+    else:
+        steps.append(("spline", method))
+    interpolator = Chain(steps)
+    interpolator.fit(coordinates, data[name])
+
+    projected = interpolator.grid(
+        region=region,
+        spacing=spacing,
+        data_names=kwargs.pop("data_names", [name]),
+        **kwargs
+    )
+    if method not in ["linear", "cubic"]:
+        projected = convexhull_mask(coordinates, grid=projected)
+    return projected[name]