diff --git a/doc/api/index.rst b/doc/api/index.rst index cce3f68f19d..6fad08b8165 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -155,6 +155,7 @@ Operations on raster data grdhisteq.equalize_grid grdhisteq.compute_bins grdlandmask + grdpaste grdproject grdsample grdtrack diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 652c6c78712..f10c4fad8a5 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -45,6 +45,7 @@ grdhisteq, grdinfo, grdlandmask, + grdpaste, grdproject, grdsample, grdtrack, diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index d91713ace78..3c3fb2c39e3 100644 --- a/pygmt/src/__init__.py +++ b/pygmt/src/__init__.py @@ -25,6 +25,7 @@ from pygmt.src.grdimage import grdimage from pygmt.src.grdinfo import grdinfo from pygmt.src.grdlandmask import grdlandmask +from pygmt.src.grdpaste import grdpaste from pygmt.src.grdproject import grdproject from pygmt.src.grdsample import grdsample from pygmt.src.grdtrack import grdtrack diff --git a/pygmt/src/grdpaste.py b/pygmt/src/grdpaste.py new file mode 100644 index 00000000000..6c9530927ce --- /dev/null +++ b/pygmt/src/grdpaste.py @@ -0,0 +1,77 @@ +""" +grdpaste - Join two grids along their common edge. +""" + +from typing import Literal + +import xarray as xr +from pygmt._typing import PathLike +from pygmt.alias import AliasSystem +from pygmt.clib import Session +from pygmt.helpers import build_arg_list, fmt_docstring, use_alias + + +@fmt_docstring +@use_alias(f="coltypes") +def grdpaste( + grid1: PathLike | xr.DataArray, + grid2: PathLike | xr.DataArray, + outgrid: PathLike | None = None, + verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"] + | bool = False, + **kwargs, +) -> xr.DataArray | None: + """ + Join two grids along their common edge. + + Requires GMT dev version (>=6.6.0). + + Combine ``grid1`` and ``grid2`` into a single grid by pasting them together along + their common edge. The two input grids must have the same grid spacings and + registration, and must have one edge in common. If in doubt, check with + :func:`pygmt.grdinfo` and use :func:`pygmt.grdcut` and/or :func:`pygmt.grdsample` if + necessary to prepare the edge joint. Note: For geographical grids, you may have to + use ``coltypes`` to handle periodic longitudes unless the input grids are properly + recognized as such via their meta-data. For stitching multiple grids, see + ``grdblend`` (not implemented in PyGMT yet) instead. + + Full GMT docs at :gmt-docs:`grdpaste.html`. + + $aliases + - G = outgrid + - V = verbose + + Parameters + ---------- + grid1 + grid2 + The two grids to be pasted. Each can be a file name or a + :class:`xarray.DataArray` object. + $outgrid + $verbose + $coltypes + + Returns + ------- + ret + 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 the file set by + ``outgrid``) + """ + aliasdict = AliasSystem().add_common(V=verbose) + aliasdict.merge(kwargs) + + with Session() as lib: + with ( + lib.virtualfile_in(check_kind="raster", data=grid1) as vingrd1, + lib.virtualfile_in(check_kind="raster", data=grid2) as vingrd2, + lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd, + ): + aliasdict["G"] = voutgrd + lib.call_module( + module="grdpaste", + args=build_arg_list(aliasdict, infile=[vingrd1, vingrd2]), + ) + return lib.virtualfile_to_raster(vfname=voutgrd, outgrid=outgrid) diff --git a/pygmt/tests/test_grdpaste.py b/pygmt/tests/test_grdpaste.py new file mode 100644 index 00000000000..ce768ef12bc --- /dev/null +++ b/pygmt/tests/test_grdpaste.py @@ -0,0 +1,121 @@ +""" +Test pygmt.grdpaste. +""" + +import pytest +import xarray as xr +from packaging.version import Version +from pygmt import grdcut, grdpaste +from pygmt.clib import __gmt_version__ +from pygmt.helpers import GMTTempFile +from pygmt.helpers.testing import load_static_earth_relief + + +@pytest.fixture(scope="module", name="grid") +def fixture_grid(): + """ + Load the grid data from the sample earth_relief file. + """ + return load_static_earth_relief() + + +@pytest.fixture(scope="module", name="grid_top") +def fixture_grid_top(grid): + """ + Load the top part of the grid data from the sample earth_relief file. + """ + return grdcut(grid, region=[-53, -49, -19, -16]) + + +@pytest.fixture(scope="module", name="grid_bottom") +def fixture_grid_bottom(grid): + """ + Load the bottom part of the grid data from the sample earth_relief file. + """ + return grdcut(grid, region=[-53, -49, -22, -19]) + + +def test_grdpaste_file_in_file_out(grid): + """ + Test grdpaste with file input and file output. + """ + with ( + GMTTempFile(suffix=".nc") as tmp1, + GMTTempFile(suffix=".nc") as tmp2, + GMTTempFile(suffix=".nc") as tmpout, + ): + grdcut(grid, region=[-53, -49, -19, -16], outgrid=tmp1.name) + grdcut(grid, region=[-53, -49, -22, -19], outgrid=tmp2.name) + result = grdpaste(grid1=tmp1.name, grid2=tmp2.name, outgrid=tmpout.name) + assert result is None # grdpaste returns None if output to a file + temp_grid = xr.load_dataarray(tmpout.name, engine="gmt", raster_kind="grid") + assert isinstance(temp_grid, xr.DataArray) + assert temp_grid.shape == (6, 4) + # Check that the result has the expected min and max values + assert temp_grid.min().values == 345.5 + assert temp_grid.max().values == 886.0 + + +# TODO(GMT>6.6.0): Remove the xfail marker for GMT<6.7. +@pytest.mark.xfail( + condition=Version(__gmt_version__) > Version("6.6.0"), + reason="Upstream bug fixed in https://github.com/GenericMappingTools/gmt/pull/8901", +) +def test_grdpaste_file_in_xarray_out(grid): + """ + Test grdpaste with file input and xarray output. + """ + with ( + GMTTempFile(suffix=".nc") as tmp1, + GMTTempFile(suffix=".nc") as tmp2, + ): + grdcut(grid, region=[-53, -49, -19, -16], outgrid=tmp1.name) + grdcut(grid, region=[-53, -49, -22, -19], outgrid=tmp2.name) + result = grdpaste(grid1=tmp1.name, grid2=tmp2.name) + assert isinstance(result, xr.DataArray) + assert result.shape == (6, 4) + # Check that the result has the expected min and max values + assert result.min().values == 345.5 + assert result.max().values == 886.0 + + +# TODO(GMT>6.6.0): Remove the xfail marker. +@pytest.mark.xfail( + condition=Version(__gmt_version__) <= Version("6.6.0"), + reason="Upstream bug fixed in https://github.com/GenericMappingTools/gmt/pull/8901", +) +def test_grdpaste(grid_top, grid_bottom): + """ + Test grdpaste by pasting two grids together along their common edge. + """ + # Paste the two grids back together + result = grdpaste(grid1=grid_top, grid2=grid_bottom) + # Check that the result is a DataArray + assert isinstance(result, xr.DataArray) + # Check that the result has the expected shape + # grid_top has 3x4, grid_bottom has 3x4, so result should have 6x4 + assert result.shape == (6, 4) + # Check that the result has the expected min and max values + assert result.min().values == 345.5 + assert result.max().values == 886.0 + + +# TODO(GMT>6.6.0): Remove the xfail marker. +@pytest.mark.xfail( + condition=Version(__gmt_version__) <= Version("6.6.0"), + reason="Upstream bug fixed in https://github.com/GenericMappingTools/gmt/pull/8901", +) +def test_grdpaste_outgrid(grid_top, grid_bottom): + """ + Test grdpaste with outgrid parameter. + """ + # Paste the two grids back together and save to file + with GMTTempFile(suffix=".nc") as tmpfile: + result = grdpaste(grid1=grid_top, grid2=grid_bottom, outgrid=tmpfile.name) + assert result is None # grdpaste returns None if output to a file + temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid") + assert isinstance(temp_grid, xr.DataArray) + assert temp_grid.shape == (6, 4) + # Check that the result has the expected min and max values + assert temp_grid.min().values == 345.5 + assert temp_grid.max().values == 886.0