Repository: fatiando/verde
Pull request #246 (2020-03-26 11:14:00)

Add function to project 2D gridded data
The `project_grid` function transforms a grid to the given projection. It re-samples the data using `ScipyGridder` (by default) and runs a blocked mean (optional) to avoid aliasing when the points aren't evenly distributed in the projected coordinates (like in polar projections). Applies a `convexhul_mask` to the grid always to avoid extrapolation to points that had no original data (this is done by scipy for linear and cubic interpolation already). 

The function only works for `xarray.DataArray`s (not `Dataset`) because I couldn't find an easy way to guarantee that a variables have the same dimensions. We can add that option later on if needed. It wouldn't break the API so we can use this as a minimum viable product for now.

Move the new function and `project_region` to the new module `verde/projections.py`.

**Reminders**:

- [x] Run `make format` and `make check` to make sure the code follows the style guide.
- [x] Add tests for new features or tests that would have caught the bug that you're fixing.
- [x] Add new public functions/methods/classes to `doc/api/index.rst` and the base `__init__.py` file for the package.
- [x] Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
- [x] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
- [x] Add your full name, affiliation, and [ORCID](https://orcid.org) (optional) to the `AUTHORS.md` file (if you haven't already) in case you'd like to be listed as an author on the [Zenodo](https://zenodo.org/communities/fatiando) archive of the next release.


leouieda: This is what the gallery example looks like:

```python
import numpy as np
import matplotlib.pyplot as plt
import pyproj
import verde as vd

# We'll use synthetic data near to 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 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()
```

![sphx_glr_project_grid_001](https://user-images.githubusercontent.com/290082/77688516-19fdd880-6f98-11ea-9904-092cfd89568b.png)

leouieda: @djhoese I was hoping to get your input on this PR as it overlaps with pyresample quite a bit (reviving our discussion in #108). This is a lot simpler than what pyresample can do but it's often what I need and not worth the effort of setting up the `AreaDefinition` etc in pyresample. But I'd like to bring this to your attention as a possible way to collaborate :slightly_smiling_face: 

I'm particularly interested in the parallel kdtree and dask support in pyresample (but that's beyond this PR).
jessepisel: @leouieda This looks like a great contribution. I will go ahead with a code review tomorrow and see what I can find.
leouieda: Thanks @jessepisel, I really appreciate it :slightly_smiling_face: 
leouieda: > I assume is because it's straightforward enough to not warrant testing

You'd be surprised! These little functions to get the coordinates right were some of the trickiest ones in the whole project. The tests for that function are the examples in the docstring (doctests). I like those if the function outputs are simple enough that they can be checked against expected results. Then it serves the dual purpose of documentation + test. 

Thanks for the review! Merging this in so it can be part of 1.4.0.

Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index bcc49d8..39a5b53 100644
--- a/doc/api/index.rst
+++ b/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
 -------
 
diff --git a/examples/project_grid.py b/examples/project_grid.py
new file mode 100644
index 0000000..5e80b49
--- /dev/null
+++ b/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()
diff --git a/verde/__init__.py b/verde/__init__.py
index a77e65c..c0b5cf5 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -12,7 +12,6 @@ from .coordinates import (
     profile_coordinates,
     get_region,
     pad_region,
-    project_region,
     longitude_continuity,
 )
 from .mask import distance_mask, convexhull_mask
@@ -26,6 +25,7 @@ from .chain import Chain
 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):
diff --git a/verde/coordinates.py b/verde/coordinates.py
index ea6da7c..b43ebd3 100644
--- a/verde/coordinates.py
+++ b/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.
diff --git a/verde/projections.py b/verde/projections.py
new file mode 100644
index 0000000..a8a580f
--- /dev/null
+++ b/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]


Test patch:
diff --git a/verde/tests/test_projections.py b/verde/tests/test_projections.py
new file mode 100644
index 0000000..ed6af18
--- /dev/null
+++ b/verde/tests/test_projections.py
@@ -0,0 +1,115 @@
+"""
+Test the projection functions.
+"""
+import numpy.testing as npt
+import numpy as np
+import xarray as xr
+import pytest
+
+from ..scipygridder import ScipyGridder
+from ..projections import project_grid
+
+
+def projection(longitude, latitude):
+    "Dummy projection"
+    return longitude ** 2, latitude ** 2
+
+
+@pytest.mark.parametrize(
+    "method",
+    ["nearest", "linear", "cubic", ScipyGridder("nearest")],
+    ids=["nearest", "linear", "cubic", "gridder"],
+)
+def test_project_grid(method):
+    "Use a simple projection to test that the output is as expected"
+    shape = (50, 40)
+    lats = np.linspace(2, 10, shape[1])
+    lons = np.linspace(-10, 2, shape[0])
+    data = np.ones(shape, dtype="float")
+    grid = xr.DataArray(data, coords=[lons, lats], dims=("latitude", "longitude"))
+    proj = project_grid(grid, projection, method=method)
+    assert proj.dims == ("northing", "easting")
+    assert proj.name == "scalars"
+    assert proj.shape == shape
+    # Check the grid spacing is constant
+    spacing_east = proj.easting[1:] - proj.easting[0:-1]
+    npt.assert_allclose(spacing_east, spacing_east[0])
+    spacing_north = proj.northing[1:] - proj.northing[0:-1]
+    npt.assert_allclose(spacing_north, spacing_north[0])
+    # Check that the values are all 1
+    npt.assert_allclose(proj.values[~np.isnan(proj.values)], 1)
+
+
+def test_project_grid_name():
+    "Check that grid name is kept"
+    shape = (50, 40)
+    lats = np.linspace(2, 10, shape[1])
+    lons = np.linspace(-10, 2, shape[0])
+    data = np.ones(shape, dtype="float")
+    grid = xr.DataArray(
+        data, coords=[lons, lats], dims=("latitude", "longitude"), name="yara"
+    )
+    proj = project_grid(grid, projection)
+    assert proj.name == "yara"
+    assert proj.dims == ("northing", "easting")
+    assert proj.shape == shape
+    # Check the grid spacing is constant
+    spacing_east = proj.easting[1:] - proj.easting[0:-1]
+    npt.assert_allclose(spacing_east, spacing_east[0])
+    spacing_north = proj.northing[1:] - proj.northing[0:-1]
+    npt.assert_allclose(spacing_north, spacing_north[0])
+    # Check that the values are all 1
+    npt.assert_allclose(proj.values[~np.isnan(proj.values)], 1)
+
+
+@pytest.mark.parametrize("antialias", [True, False])
+def test_project_grid_antialias(antialias):
+    "Check if antialias is being used"
+    shape = (50, 40)
+    lats = np.linspace(2, 10, shape[1])
+    lons = np.linspace(-10, 2, shape[0])
+    data = np.ones(shape, dtype="float")
+    grid = xr.DataArray(data, coords=[lons, lats], dims=("latitude", "longitude"))
+    proj = project_grid(grid, projection, antialias=antialias)
+    if antialias:
+        assert "BlockReduce" in proj.attrs["metadata"]
+    else:
+        assert "BlockReduce" not in proj.attrs["metadata"]
+    assert proj.dims == ("northing", "easting")
+    assert proj.name == "scalars"
+    assert proj.shape == shape
+    # Check the grid spacing is constant
+    spacing_east = proj.easting[1:] - proj.easting[0:-1]
+    npt.assert_allclose(spacing_east, spacing_east[0])
+    spacing_north = proj.northing[1:] - proj.northing[0:-1]
+    npt.assert_allclose(spacing_north, spacing_north[0])
+    # Check that the values are all 1
+    npt.assert_allclose(proj.values[~np.isnan(proj.values)], 1)
+
+
+def test_project_grid_fails_dataset():
+    "Should raise an exception when given a Datatset"
+    shape = (50, 40)
+    lats = np.linspace(2, 10, shape[1])
+    lons = np.linspace(-10, 2, shape[0])
+    data = np.ones(shape, dtype="float")
+    grid = xr.DataArray(data, coords=[lons, lats], dims=("latitude", "longitude"))
+    grid = grid.to_dataset(name="scalars")
+    with pytest.raises(ValueError):
+        project_grid(grid, projection)
+
+
+@pytest.mark.parametrize("ndims", [1, 3])
+def test_project_grid_fails_dimensions(ndims):
+    "Should raise an exception when given more or less than 2 dimensions"
+    shape = (10, 20, 12)
+    coords = [
+        np.linspace(-10, 2, shape[0]),
+        np.linspace(2, 10, shape[1]),
+        np.linspace(0, 100, shape[2]),
+    ]
+    dims = ("height", "latitude", "longitude")
+    data = np.ones(shape[:ndims], dtype="float")
+    grid = xr.DataArray(data, coords=coords[:ndims], dims=dims[:ndims])
+    with pytest.raises(ValueError):
+        project_grid(grid, projection)