grdtrack/surface/others?, incompatible data types

[ashtonflinders]

Although the docs says pygmt.surface returns a xarray.DataArray, it actually returns an xarray.Dataset. This is incompatible to what grdtrack expects as input, making it not possible to sample an internally created grid unless you write it out to file, then reread it, or do some finagling.

import pygmt
bathy= pygmt.datasets.load_sample_bathymetry()

xmin = bathy['longitude'].min()
xmax = bathy['longitude'].max()
ymin = bathy['latitude'].min()
ymax = bathy['latitude'].max()

R = '/'.join([str(xmin),str(xmax),str(ymin),str(ymax)])
I='.02'

x = bathy['longitude']
y = bathy['latitude']
z = bathy['bathymetry']

DEPTH_SURF = pygmt.surface(x = x, y = y, z = z, R = R, I = I) 

type(DEPTH_SURF)
xarray.core.dataset.Dataset

pygmt.grdtrack(bathy,DEPTH_SURF,newcolname='sampled')

GMTInvalidInput: Unrecognized data type xarray.core.dataset.Dataset

The work around is to convert the dataset to a dataarray;

depth_array=DEPTH_SURF.to_array(dim='z')

Although this doesn't work directly, since it incorrectly broadcasts the array size;

pygmt.grdtrack(bathy,depth_array,newcolname='sampled')

GMTInvalidInput: Invalid number of grid dimensions '3'. Must be 2.

depth_array.values.shape
(1, 501, 486)

To get it to work, you have to use a slice of the dataarray;

pygmt.grdtrack(bathy,depth_array[0],newcolname='sampled')

       longitude  latitude  bathymetry      sampled
0      245.00891  27.49555        -636  -690.523706
1      245.01201  27.49286        -655  -697.038135
2      245.01512  27.49016        -710  -703.130958
...          ...       ...         ...          ...
82968  245.00681  24.61633       -3852 -3861.402543
82969  245.00337  24.61995       -3889 -3869.346273

Unless I'm missing something obvious here and this is the expected behavior?

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

[weiji14]

Oops, you're right, it should return an xr.DataArray instead. I'll work on a fix for it now.

In the meantime, you can use DEPTH_SURF.z to access the xr.DataArray inside the xr.Dataset.