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..6fa95ef564a --- /dev/null +++ b/pygmt/src/grdpaste.py @@ -0,0 +1,94 @@ +""" +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 + +__doctest_skip__ = ["grdpaste"] + + +@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``) + + Example + ------- + >>> import pygmt + >>> # Create two grids with a common edge + >>> # Grid 1: longitude range of 10° E to 20° E, latitude range of 15° N to 25° N + >>> grid1 = pygmt.datasets.load_earth_relief( + ... resolution="30m", region=[10, 20, 15, 25] + ... ) + >>> # Grid 2: longitude range of 10° E to 20° E, latitude range of 10° N to 15° N + >>> grid2 = pygmt.datasets.load_earth_relief( + ... resolution="30m", region=[10, 20, 10, 15] + ... ) + >>> # Paste the two grids together along their common edge (15° N) + >>> new_grid = pygmt.grdpaste(grid1=grid1, grid2=grid2) + """ + 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..9e6e4170708 --- /dev/null +++ b/pygmt/tests/test_grdpaste.py @@ -0,0 +1,119 @@ +""" +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, grid_top, grid_bottom): + """ + 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 == min( + grid_top.min().values, grid_bottom.min().values + ) + assert temp_grid.max().values == max( + grid_top.max().values, grid_bottom.max().values + ) + + +# TODO(GMT>=6.7): Remove the xfail marker for GMT<6.7. +@pytest.mark.xfail( + condition=Version(__gmt_version__) < Version("6.7"), + reason="Requires GMT dev version (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) + + +# TODO(GMT>=6.7): Remove the xfail marker for GMT<6.7. +@pytest.mark.xfail( + condition=Version(__gmt_version__) < Version("6.7"), + reason="Requires GMT dev version (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 == min(grid_top.min().values, grid_bottom.min().values) + assert result.max().values == max(grid_top.max().values, grid_bottom.max().values) + + +# TODO(GMT>=6.7): Remove the xfail marker for GMT<6.7. +@pytest.mark.xfail( + condition=Version(__gmt_version__) < Version("6.7"), + reason="Requires GMT dev version (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)