pygmt.grdpaste: Join two grids along their common edge

[Chuan1937]

pygmt.grdpaste: Join two grids along their common edge

Preview: https://pygmt-dev--4399.org.readthedocs.build/en/4399/api/generated/pygmt.grdpaste.html#pygmt.grdpaste

[seisman]

Also need to add it to doc/api/index.rst.

[Chuan1937]

Could you also please check the min and max values in the result?
Must use temporary files for xarray.DataArray inputs because virtualfile_in returns all zeros when handling grdpaste.

[seisman]

Could you also please check the min and max values in the result? Must use temporary files for xarray.DataArray inputs because virtualfile_in returns all zeros when handling grdpaste.

I think this points to an upstream GMT bug, which should be fixed.

[seisman]

Could you please add one more test with grid1 as file and grid2 as xr.dataarray？

[Chuan1937]

Could you please add one more test with grid1 as file and grid2 as xr.dataarray？

When the first input is a real file path and the second is a virtual file (from a DataArray), a GMT_GRID_READ_ERROR is raised. However, if the first input is a virtual file and the second is a real file, it works normally.

success:

def test_grdpaste_grid1_file_grid2_xarray(grid, grid_bottom):
    """
    Test grdpaste with grid1 as file input and grid2 as xarray input.
    """
    with GMTTempFile(suffix=".nc") as tmp1:
        grdcut(grid, region=[-53, -49, -19, -16], outgrid=tmp1.name)
        result = grdpaste(grid1=tmp1.name, grid2=grid_bottom)
        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

failure:

def test_grdpaste_grid1_xarray_grid2_file(grid, grid_top):
    """
    Test grdpaste with grid1 as xarray input and grid2 as file input.
    """
    with GMTTempFile(suffix=".nc") as tmp2:
        grdcut(grid, region=[-53, -49, -22, -19], outgrid=tmp2.name)
        result = grdpaste(grid1=grid_top, 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

A temporary workaround is to detect mixed inputs in PyGMT and load the files as DataArrays, making both inputs virtual files:

 if {kind1, kind2} == {"file", "grid"}:
            if kind1 == "file":
                grid1 = xr.load_dataarray(grid1, engine="gmt", raster_kind="grid")
            else:
                grid2 = xr.load_dataarray(grid2, engine="gmt", raster_kind="grid")

This may be an issue with upstream GMT. Should we use the workaround or report the problem to upstream?
@seisman

[seisman]

@Chuan1937 I've added a section "Notes on the input/output" in #4395. Could you please complete the summary table so that we can better understand the current behavior/limitations?

[Chuan1937]

@seisman
I've completed "Notes on the input/output".

[seisman]

Thanks. Could you please add the tests in #4399 (comment)? You can mark the failing one as xfail.

[seisman]

Thanks. Could you please add the tests in #4399 (comment)? You can mark the failing one as xfail.

Wait. Maybe we should make it simple by only allowing grid1 and grid2 to be the same type (file or xarray.DataArray).

[Chuan1937]

Thanks. Could you please add the tests in #4399 (comment)? You can mark the failing one as xfail.

Wait. Maybe we should make it simple by only allowing grid1 and grid2 to be the same type (file or xarray.DataArray).

I agree. It seems there are still quite a few issues with GMT's external interfaces lately. I'll add this test first.