grdtrack produces incorrect results for xarray.DataArray grid

[seisman]

Description of the problem

grdtrack works well for netCDF file, but produces random results for xarray.DataArray grid. If I understand #1304 correctly, it works on Linux, but not on macOS and Windows. I believe it's an upstream GMT bug.

Full code that generated the error

from pygmt import grdtrack, grdcut
from pygmt.datasets import load_earth_relief

# load the earth_relief_01d_g file
grid = load_earth_relief(registration="gridline")

# prepare a subgrid in xarray.DataArray
dataarray = grid.sel(lat=slice(-49, -42), lon=slice(-118, -107))
# prepare a subgrid in netcdf
grdcut(grid, region=[-118, -107, -49, -42], outgrid="earth_subset.nc")

# run grdtrack on the netCDF file
output = grdtrack(points="@ridge.txt", grid="earth_subset.nc", newcolname="bathymetry")
print(output)

# run grdtrack on the xarray.DataArray
output = grdtrack(points="@ridge.txt", grid=dataarray, newcolname="bathymetry")
print(output)

Full output

The first output is correct, and the second output is wrong.

grdtrack [WARNING]: Some input points were outside the grid domain(s).
          0        1            2
0 -110.9536 -42.2489 -2797.394987
1 -111.3500 -42.2033 -2715.359945
2 -111.7309 -43.1706 -2711.704748
3 -112.2196 -44.3819 -2787.948153
4 -112.5505 -45.4878 -2746.303734
5 -112.6736 -45.8347 -2685.852247
6 -112.4235 -45.8561 -2785.099618
7 -112.9507 -47.2270 -2672.932596
8 -113.3124 -48.3208 -2706.345811
grdtrack [WARNING]: Some input points were outside the grid domain(s).
          0        1             2
0 -110.9536 -42.2489 -2.974656e+03
1 -111.3500 -42.2033 -2.859961e+03
2 -111.7309 -43.1706 -2.711705e+03
3 -112.2196 -44.3819 -2.787948e+03
4 -112.5505 -45.4878 -2.746304e+03
5 -112.6736 -45.8347 -2.685852e+03
6 -112.4235 -45.8561 -2.785100e+03
7 -112.9507 -47.2270 -2.672933e+03
8 -113.3124 -48.3208  1.123369e+28

System information

Please paste the output of python -c "import pygmt; pygmt.show_versions()":

PyGMT information:
  version: v0.3.2.dev168+ga8472925
System information:
  python: 3.9.2 | packaged by conda-forge | (default, Feb 21 2021, 05:02:20)  [Clang 11.0.1 ]
  executable: /Users/seisman/.miniconda/envs/pygmt/bin/python
  machine: macOS-11.2.3-x86_64-i386-64bit
Dependency information:
  numpy: 1.20.1
  pandas: 1.2.3
  xarray: 0.18.0
  netCDF4: 1.5.6
  packaging: 20.9
  ghostscript: 9.54.0
  gmt: 6.2.0_922af84_2021.03.17
GMT library information:
  binary dir: /Users/seisman/.miniconda/envs/pygmt/bin
  cores: 8
  grid layout: rows
  library path: /Users/seisman/opt/GMT-master/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/seisman/opt/GMT-master/lib/gmt/plugins
  share dir: /Users/seisman/opt/GMT-master/share
  version: 6.2.0

@meghanrjones @PaulWessel Please see if you can debug it.

[PaulWessel]

I get the first output but on the 2nd one I get a crash:

>>> output = grdtrack(points="@ridge.txt", grid=dataarray, newcolname="bathymetry")
Assertion failed: (node < G->header->size), function gmt_bcr_get_z_fast, file /Users/pwessel/GMTdev/gmt-dev/src/gmt_bcr.c, line 256.
Abort trap: 6

Key problem is that grdtrack assumes bicubic sampling which requires padding and your matrix has none. The other place this is an issue is in grdimage and there we create a padded grid if resampling is needed. Not sure if we should do the same here, but if we don't then you will never get to use bicubic. The alternative is to drop to bilinear and then no pads are needed. Perhaps the right thing to do is to see if -n is set of defaults to anything but bilinear and you have no pad then we must create that padded grid.

[PaulWessel]

I also wonder if the more general solution to these are to possibly add a flag to the mode passed to GMT_Read_Data that means "this grid needs to have 2 boundary pads and if not then make it so". This flag would not do much for the case of reading from a data file, but in the case of a read-only input matrix with no pad it would basically allocate the grid, copy in the matrix, and from there on out we are working with the GMT-allocated grid. So, this flag would be set based on need, as follows:

1. grdgrack: If -n implies the need for 2 rows/cols, pass GMT_GRID_NEEDS_PAD else not.
2. grdimage: If need_to_project = (gmt_M_is_nonlinear_graticule (GMT) || Ctrl->E.dpi > 0); is true then pass GMT_GRID_NEEDS_PAD in any GMT_Read_Data calls, else not.

I think this is the best solution rather than adding ad-hoc stuff like currently done in gmt_grd_project and what would be needed in grdtrack, etc. However, this may not be a super-quick fix, so perhaps in the interest of 6.2.0 we could temporarily impose -nl if no pad on the incoming grid? If so, should this be done in PyGMT?

[PaulWessel]

If you have a way of passing -nl to grdtrack in pygmt then see if the script runs as is.

[seisman]

Yes, adding n="l" gives the same results for the two grdtrack calls.

I also wonder if the more general solution to these are to possibly add a flag to the mode passed to GMT_Read_Data that means "this grid needs to have 2 boundary pads and if not then make it so".

Does it means that wrappers like PyGMT must know the technical details to decide if GMT_GRID_NEEDS_PADS is necessary?

[PaulWessel]

No, only the C modules need to know what they need, so easy for the wrappers: No changes (after implementing the GMT_GRID_NEEDS_PADS). I think @joa-quim may agree that this is a good enhancement. So in that case the GMT_GRID_NEEDS_PADS is not an API addition but just another internal developer constant.

[seisman]

No, only the C modules need to know what they need, so easy for the wrappers: No changes (after implementing the GMT_GRID_NEEDS_PADS).

It sounds good.

[PaulWessel]

So I think until then the simple thing is to check in grdtrack that if pad is < 2 then we set the -n interpolant to bilinear.

[seisman]

FYI, the issue has been temporarily fixed in GenericMappingTools/gmt#5283 and should be available in GMT 6.2.0.