Passing xarray.DataArray to grdfilter gives incorrect results

[seisman]

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.

[seisman]

@PaulWessel It looks like a GMT bug. Please see if you can fix it when you have time.

[PaulWessel]

To keep it as simple as possible, can you demonstrate the same mismatch with -Fg300 or even better, -Fb300 ?

[seisman]

CLI:

gmt grdfilter @earth_relief_01d_p -Fb300 -D4 -Goutput.nc
gmt grdinfo -C output.nc
output.nc	-180	180	-90	90	-6470.13720703	5266.73339844	1	1	360180	1	1

Python:

from pygmt import grdfilter
from pygmt.datasets import load_earth_relief

# passing a grid file
result1 = grdfilter("@earth_relief_01d_p", filter="b300", 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="b300", distance="4")
print("result2:", result2.data.min(), result2.data.max())

the output is even worse:

result1: -6470.137 5266.7334
result2: nan nan

[PaulWessel]

How would you do the equivalent of

gmt grdcut @earth_relief_01d_p -R0/10/0/10 -Gt.nc
gmt grdfilter t.grd -Fb300 -D4 -Goutput.nc

in python, and does it differ then too (trying to make as small test as possible)

[seisman]

import pygmt

grid = pygmt.grdcut("@earth_relief_01d_p", region=[0, 10, 0, 10])
print(grid.data.max(), grid.data.min())

result = pygmt.grdfilter(grid, filter='b300', distance='4')
print(result.data.max(), result.data.min())

output:

821.5 -4852.0
nan nan

[seisman]

Ping @PaulWessel just in case you forget it.

[PaulWessel]

Was starting to look using my new M1 laptop but first have to install miniconda and that got me this far:

bash ~/Downloads/Miniconda3-py39_4.9.2-MacOSX-x86_64.sh
WARNING:
Your operating system appears not to be 64-bit, but you are trying to
install a 64-bit version of Miniconda3.
Are sure you want to continue the installation? [yes|no]

Some googling talks about condo running in emulation mode but not sure about miniconda. I guess this is all work in progress. Not sure how to proceed. Perhaps just not do this on laptop but do it on the old Mac Pro Intel...

[PaulWessel]

Hi @seisman, is this still failing after the last API fixes today/yesterday?

[seisman]

It's still failing.

[PaulWessel]

Odd, I just stepped through the whole grdfilter example above (the 10x10 case) and I got the same result as CLI. Sure you rebuilt etc?

[seisman]

I think I'm using the latest GMT version (6.2.0_c673720_2021.04.05) and all three examples still fail.

What happens if you just run the examples, not in Xcode debugging.

[PaulWessel]

(base) pwessel@macnut:~/GMTdev/pygmt-> python
Python 3.8.3 (default, May 19 2020, 13:54:14) 
[Clang 10.0.0 ] :: Anaconda, Inc. on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pygmt
>>> 
>>> grid = pygmt.grdcut("@earth_relief_01d_p", region=[0, 10, 0, 10])
>>> print(grid.data.max(), grid.data.min())
821.5 -4852.0
>>> 
>>> result = pygmt.grdfilter(grid, filter='b300', distance='4')
>>> print(result.data.max(), result.data.min())
599.0929 -4787.664

Same result with release build or xbuild.