grdtrack/surface/others?, incompatible data types #407
Description
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?