Wrap grdproject (GenericMappingTools#1377)

Co-authored-by: Wei Ji <[email protected]>
Co-authored-by: Dongdong Tian <[email protected]>

doc/api/index.rst

@@ -94,6 +94,7 @@ Operations on grids:
     grdfilter
     grdgradient
     grdlandmask
+    grdproject
     grdsample
     grdtrack
     xyz2grd

pygmt/__init__.py

@@ -40,6 +40,7 @@
     grdgradient,
     grdinfo,
     grdlandmask,
+    grdproject,
     grdsample,
     grdtrack,
     info,

pygmt/src/__init__.py

@@ -19,6 +19,7 @@
 from pygmt.src.grdimage import grdimage
 from pygmt.src.grdinfo import grdinfo
 from pygmt.src.grdlandmask import grdlandmask
+from pygmt.src.grdproject import grdproject
 from pygmt.src.grdsample import grdsample
 from pygmt.src.grdtrack import grdtrack
 from pygmt.src.grdview import grdview

pygmt/src/grdproject.py

@@ -0,0 +1,93 @@
+"""
+grdproject - Forward and inverse map transformation of grids.
+"""
+
+import xarray as xr
+from pygmt.clib import Session
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+
+
+@fmt_docstring
+@use_alias(
+    G="outgrid",
+    J="projection",
+    I="inverse",
+    R="region",
+    V="verbose",
+    n="interpolation",
+    r="registration",
+)
+@kwargs_to_strings(R="sequence")
+def grdproject(grid, **kwargs):
+    r"""
+    Change projection of gridded data between geographical and rectangular.
+
+    This module will project a geographical gridded data set onto a
+    rectangular grid. If ``inverse`` is ``True``, it will project a
+    rectangular coordinate system to a geographic system. To obtain the value
+    at each new node, its location is inversely projected back onto the input
+    grid after which a value is interpolated between the surrounding input
+    grid values. By default bi-cubic interpolation is used. Aliasing is
+    avoided by also forward projecting the input grid nodes. If two or more
+    nodes are projected onto the same new node, their average will dominate in
+    the calculation of the new node value. Interpolation and aliasing is
+    controlled with the ``interpolation`` option. The new node spacing may be
+    determined in one of several ways by specifying the grid spacing, number
+    of nodes, or resolution. Nodes not constrained by input data are set to
+    NaN. The ``region`` parameter can be used to select a map region larger or
+    smaller than that implied by the extent of the grid file.
+
+    {aliases}
+
+    Parameters
+    ----------
+    grid : str or xarray.DataArray
+        The file name of the input grid or the grid loaded as a DataArray.
+    outgrid : str or None
+        The name of the output netCDF file with extension .nc to store the grid
+        in.
+    inverse : bool
+        When set to ``True`` transforms grid from rectangular to
+        geographical [Default is False].
+    {J}
+    {R}
+    {V}
+    {n}
+    {r}
+
+    Returns
+    -------
+    ret: xarray.DataArray or None
+        Return type depends on whether the ``outgrid`` parameter is set:
+
+        - :class:`xarray.DataArray` if ``outgrid`` is not set
+        - None if ``outgrid`` is set (grid output will be stored in file set by
+          ``outgrid``)
+    """
+    if "J" not in kwargs.keys():
+        raise GMTInvalidInput("The projection must be specified.")
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        with Session() as lib:
+            file_context = lib.virtualfile_from_data(check_kind="raster", data=grid)
+            with file_context as infile:
+                if "G" not in kwargs.keys():  # if outgrid is unset, output to tempfile
+                    kwargs.update({"G": tmpfile.name})
+                outgrid = kwargs["G"]
+                arg_str = " ".join([infile, build_arg_string(kwargs)])
+                lib.call_module("grdproject", arg_str)
+
+        if outgrid == tmpfile.name:  # if user did not set outgrid, return DataArray
+            with xr.open_dataarray(outgrid) as dataarray:
+                result = dataarray.load()
+                _ = result.gmt  # load GMTDataArray accessor information
+        else:
+            result = None  # if user sets an outgrid, return None
+
+        return result

pygmt/tests/test_grdproject.py

@@ -0,0 +1,59 @@
+"""
+Tests for grdproject.
+"""
+import os
+
+import numpy.testing as npt
+import pytest
+from pygmt import grdinfo, grdproject
+from pygmt.datasets import load_earth_relief
+from pygmt.exceptions import GMTInvalidInput
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="grid")
+def fixture_grid():
+    """
+    Load the grid data from the sample earth_relief file.
+    """
+    return load_earth_relief(resolution="01d", region=[-5, 5, -5, 5])
+
+
+def test_grdproject_file_out(grid):
+    """
+    grdproject with an outgrid set.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = grdproject(grid=grid, projection="M10c", outgrid=tmpfile.name)
+        assert result is None  # return value is None
+        assert os.path.exists(path=tmpfile.name)  # check that outgrid exists
+        result = grdinfo(tmpfile.name, per_column=True).strip().split()
+        npt.assert_allclose(float(result[0]), 0)  # x min
+        npt.assert_allclose(float(result[1]), 10)  # x max
+        npt.assert_allclose(float(result[2]), 0, atol=1.0e-10)  # y min
+        npt.assert_allclose(float(result[3]), 9.94585661273)  # y max
+        npt.assert_allclose(float(result[4]), -5130.48193359)  # min
+        npt.assert_allclose(float(result[5]), -152.585281372)  # max
+
+
+def test_grdproject_no_outgrid(grid):
+    """
+    Test grdproject with no set outgrid.
+    """
+    assert grid.gmt.gtype == 1  # Geographic grid
+    temp_grid = grdproject(grid=grid, projection="M10c")
+    assert temp_grid.dims == ("y", "x")
+    assert temp_grid.gmt.gtype == 0  # Rectangular grid
+    assert temp_grid.gmt.registration == 1  # Pixel registration
+    npt.assert_allclose(temp_grid.min(), -5130.482)
+    npt.assert_allclose(temp_grid.max(), -152.58528)
+    npt.assert_allclose(temp_grid.median(), -4578.288)
+    npt.assert_allclose(temp_grid.mean(), -4354.3296)
+
+
+def test_grdproject_fails(grid):
+    """
+    Check that grdproject fails correctly.
+    """
+    with pytest.raises(GMTInvalidInput):
+        grdproject(grid=grid)