pygmt: Passing xarray.DataArray to grdfilter gives incorrect results

Description of the problem

Here is the result from GMT CLI. From the output, we can know the min and max of the filtered grid, -6147.49072266 and 5164.06005859.

$ gmt grdfilter @earth_relief_01d_p -Fg600 -D4 -Goutput.nc
$ gmt grdinfo -C output.nc
output.nc	-180	180	-90	90	-6147.49072266	5164.06005859	1	1	360	180	1	1

This is how we do the same thing in Python:

from pygmt import grdfilter
from pygmt.datasets import load_earth_relief

# passing a grid file
result1 = grdfilter("@earth_relief_01d_p", filter="g600", distance="4")
print("result1:", result1.data.min(), result1.data.max())

# passing a xarray.DataArray
grid = load_earth_relief(resolution="01d", registration="pixel")
result2 = grdfilter(grid=grid, filter="g600", distance="4")
print("result2:", result2.data.min(), result2.data.max())

The output is:

result1: -6147.4907 5164.06
result2: -6147.4727 5164.1157

Obviously, passing xarray.DataArray to grdfilter gives incorrect results.

It’s still unclear to me if this is a PyGMT bug or an upstream GMT bug.

Note: the tests in test_grdfilter.py is testing the “wrong” results. We need to update the tests after figuring out the reason.