grdtrack not working with "profile" option

[19giovi87]

Hi,

I'm trying to extract the topography from a grid along a certain path, going from A to B.
For this purpose, here is the code:

region = [14.1, 15.2, 41.0, 41.6]
grid = pygmt.datasets.load_earth_relief(resolution="01s", region=region)
profile = '14.33/41.43/14.93/41.14'         # these are the coordinates of A and B

track   = pygmt.grdtrack(points = None, grid = grid,
                                          profile = profile)

This is the error I get. It occurs with any grid provided, as long as the "profile" option is used and "points" is set to None.

File "/Users/giovanni_work/anaconda3/envs/pygmt/lib/python3.9/site-packages/pygmt/helpers/utils.py", line 69, in data_kind
    raise GMTInvalidInput("No input data provided.")

GMTInvalidInput: No input data provided.

I don't know if this is a bug or I'm just wrongly using the profile/points option is grdtrack.
Any clues on how to solve this?

Thank you!

System information

PyGMT information:
  version: v0.5.0
System information:
  python: 3.9.10 | packaged by conda-forge | (main, Feb  1 2022, 21:28:27)  [Clang 11.1.0 ]
  executable: /Users/rex/anaconda3/envs/pygmt/bin/python
  machine: macOS-11.6.1-x86_64-i386-64bit
Dependency information:
  numpy: 1.22.3
  pandas: 1.4.1
  xarray: 2022.3.0
  netCDF4: 1.5.8
  packaging: 21.3
  ghostscript: 9.54.0
  gmt: 6.3.0
GMT library information:
  binary dir: /Users/giovanni_work/anaconda3/envs/pygmt/bin
  cores: 16
  grid layout: rows
  library path: /Users/giovanni_work/anaconda3/envs/pygmt/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/giovanni_work/anaconda3/envs/pygmt/lib/gmt/plugins
  share dir: /Users/giovanni_work/anaconda3/envs/pygmt/share/gmt
  version: 6.3.0

[welcome]

👋 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]

Ah yes, I think I know what's happening, this might be a bug introduced in #1189.

In the meantime, you can try using this as a workaround. You can adjust the num=250 to the number of points you want along the line:

import pandas as pd
import pygmt
import numpy as np

region = [14.1, 15.2, 41.0, 41.6]
grid = pygmt.datasets.load_earth_relief(resolution="01s", region=region)

profile = "14.33/41.43/14.93/41.14"  # these are the coordinates of A and B

points = pd.DataFrame(
    data=np.linspace(start=(14.33, 41.43), stop=(14.93, 41.14), num=250),
    columns=["x", "y"],
)

track = pygmt.grdtrack(points=points, grid=grid, newcolname="elevation")
print(track)

produces

	x	y	elevation
0	14.330000	41.430000	1294.000000
1	14.332410	41.428835	1389.962382
2	14.334819	41.427671	1445.901806
3	14.337229	41.426506	1397.092802
4	14.339639	41.425341	1472.605211
...	...	...	...
245	14.920361	41.144659	235.882504
246	14.922771	41.143494	215.088483
247	14.925181	41.142329	170.312853
248	14.927590	41.141165	158.030026
249	14.930000	41.140000	173.000000

250 rows × 3 columns

based on this chunk of code: https://github.com/weiji14/pyissm/blob/2a30cfa461b81493f41021bac0c3621a3e30c2ec/plot_figures.py#L105-L141

[19giovi87]

Ok, I see.
For now I'm going to use the workaround you proposed. Thanks!

[seisman]

Reopening the issue so that we can better track the bug report and fix it later.

[seisman]

There are two different ways to specify input track coordinates to grdtrack:

#. Pass a file or array to points
#. Set profile (grdtrack's -E option)

PR #1867 tries to fix the issue by detecting if points or profile is given. Because points is a required parameter, users have to use

pygmt.grdtrack(points=None, grid=grid, profile=profile)

which is not friendly.

I think we should reorder the function definition to grdtrack(grid, points=None, **kwargs) so that points is optional. The change should be done in a backward-compatible way and goes into the deprecatation cycle. That means #1867 can't be merged in the v0.6.1 version.

[weiji14]

Ok, considering that there is a workaround, I guess we could merge #1867 in v0.7.0. But shouldn't we start to raise a FutureWarning in v0.6.1 that the grid parameter should be first and points is second? E.g. use a modified version of the @check_data_input_order decorator used in #1479?

Or actually, can't we just change the function definition from def grdtrack(points, grid, ..., **kwargs): to def grdtrack(points=None, grid=None, ..., **kwargs):?

[seisman]

Or actually, can't we just change the function definition from def grdtrack(points, grid, ..., **kwargs): to def grdtrack(points=None, grid=None, ..., **kwargs):?

grid is a required parameter for grdtrack, and most other grd* functions have a definition like def grdxxx(grid, ..., **kwargs). I think grdtrack should also follow the same style.

But shouldn't we start to raise a FutureWarning in v0.6.1 that the grid parameter should be first and points is second?

I'm not sure if I understand you correctly. Do you mean changing the function definition or not? If we change the function definition, then the change should be in v0.7.0, not v0.6.1 following the semantic versioning. If we don't change the function definition, the FutureWarning doesn't make any sense.