Add function to project 2D gridded data

[leouieda]

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.DataArrays (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:

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• Add your full name, affiliation, and ORCID (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 archive of the next release.

[leouieda]

This is what the gallery example looks like:

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()

[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 all that I need. To project a simple 2D grid, it's often not worth the effort of setting up the AreaDefinition etc in pyresample. So I think it's worth having this utility in Verde. But I'd like to bring this to your attention as a possible way to collaborate 🙂

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 🙂

[jessepisel]

This is a great description of how it works and how the function deals with the original data. I dig it!

[leouieda]

Thanks! I did mostly because I had to dig into the GMT code in order to see what they did. It's fancier than this but Verde is much easier to install :)

[jessepisel]

This is a useful note with link to pyresample

[jessepisel]

The only thing that I could see differently is adding a plot of the data for creating the geographic grid earlier on in the example. But the more I think about it, the more I like it down here where I can visually compare the geographic and projected coordinates. I think it looks great the way it is.

[leouieda]

That would be better in a tutorial. But since the plots in the gallery appear all at the top, it wouldn't make much difference. It took me a while to think of a way to show this function in action (we don't have any grids in Verde as sample data).

[jessepisel]

This is great to have antialiasing in here

[leouieda]

I agree 🙂

[jessepisel]

Data checks look good and print useful and meaningful error messages

[jessepisel]

I think that all your changes look good. There are no issues with the builds and all the tests have passed. The only thing I noticed was there was no test for shape_to_spacing which I assume is because it's straightforward enough to not warrant testing. Only suggestions are minor. I think this is a really useful set of functions for projecting gridded data.

[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.