Segmentation fault if a grid has zero increments

[yohaimagen]

grdimage throwing a segmentation fault insted of some proper error when mistakenly swapping xarray y and x coordination.

Full code that generated the error

X, Y = np.meshgrid(np.arange(10), np.arange(10))
Z = X ** 2 + Y ** 2
da = xr.DataArray(
    data=Z,
    dims=['y', 'x'],
    coords=dict(
        x=(['x'], Y[0, :]),
        y=(['y'], X[:, 0])
    )
)

fig = pygmt.Figure()
fig.basemap(
    frame='a',
    projection=f"X10",
    region=[0, 10, 0, 10]
)
fig.grdimage(da)
fig.show()

Full error message

Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

System information

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

PyGMT information:
  version: v0.4.1
System information:
  python: 3.7.10 | packaged by conda-forge | (default, Feb 19 2021, 15:59:12)  [Clang 11.0.1 ]
  executable: /Users/yohai/anaconda3/envs/faultSlip/bin/python
  machine: Darwin-18.7.0-x86_64-i386-64bit
Dependency information:
  numpy: 1.19.1
  pandas: 1.1.1
  xarray: 0.18.2
  netCDF4: 1.5.7
  packaging: 21.0
  ghostscript: 9.50
  gmt: 6.2.0
GMT library information:
  binary dir: /Users/yohai/anaconda3/envs/faultSlip/bin
  cores: 12
  grid layout: rows
  library path: /Users/yohai/anaconda3/envs/faultSlip/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/yohai/anaconda3/envs/faultSlip/lib/gmt/plugins
  share dir: /Users/yohai/anaconda3/envs/faultSlip/share/gmt
  version: 6.2.0

[seisman]

There is something wrong with your x and y coords.

Printing the da variable gives me:

<xarray.DataArray (y: 10, x: 10)>
array([[  0,   1,   4,   9,  16,  25,  36,  49,  64,  81],
       [  1,   2,   5,  10,  17,  26,  37,  50,  65,  82],
       [  4,   5,   8,  13,  20,  29,  40,  53,  68,  85],
       [  9,  10,  13,  18,  25,  34,  45,  58,  73,  90],
       [ 16,  17,  20,  25,  32,  41,  52,  65,  80,  97],
       [ 25,  26,  29,  34,  41,  50,  61,  74,  89, 106],
       [ 36,  37,  40,  45,  52,  61,  72,  85, 100, 117],
       [ 49,  50,  53,  58,  65,  74,  85,  98, 113, 130],
       [ 64,  65,  68,  73,  80,  89, 100, 113, 128, 145],
       [ 81,  82,  85,  90,  97, 106, 117, 130, 145, 162]])
Coordinates:
  * x        (x) int64 0 0 0 0 0 0 0 0 0 0
  * y        (y) int64 0 0 0 0 0 0 0 0 0 0

as you can see, x and y coords are all zeros.

[yohaimagen]

That's true, I didn't notice. It's because of the way I'm building the DataArray with the meshgrid. Still, I'll consider this a bug because of the segmentation fault, some error should be raise.

[seisman]

I'll consider this a bug because of the segmentation fault, some error should be raise.

I agree, but I'm not sure if the error should be done by PyGMT or by GMT.

Ping @meghanrjones

[maxrjones]

I'll consider this a bug because of the segmentation fault, some error should be raise.

I agree, but I'm not sure if the error should be done by PyGMT or by GMT.

Ping @meghanrjones

I think ideally it should be caught by both so that PyGMT users have a more readable call stack and other wrappers do not get the same problem. If the same grid is a actual file rather than a remote file, GMT catches the read error properly which suggests it is an API bug.

We could pretty easily add a check for this case in the code section below. I cannot find this documented, but I think that GMT expects monotonically increasing, evenly spaced grids even though COARDS compliant grids can be monotonically increasing or decreasing. If this is the case, we could raise GMTInvalidInput if coord_inc <=0.

pygmt/pygmt/clib/conversion.py

Lines 93 to 102 in d01da5a

	for dim in grid.dims[::-1]:
	coord = grid.coords[dim].values
	coord_incs = coord[1:] - coord[0:-1]
	coord_inc = coord_incs[0]
	if not np.allclose(coord_incs, coord_inc):
	raise GMTInvalidInput(
	"Grid appears to have irregular spacing in the '{}' dimension.".format(
	dim
	)
	)

[yohaimagen]

I'll add the error.

[yohaimagen]

actually, GMT is handling well monotonically decreasing grid when passing with DataArray

see example

X, Y = np.meshgrid(np.linspace(0, -10, 10), np.linspace(0, -10, 10))
Z = X ** 2 + Y ** 2
da = xr.DataArray(
    data=Z,
    dims=['y', 'x'],
    coords=dict(
        x=(['x'], X[0, :]),
        y=(['y'], Y[:, 0])
    )
)

fig = pygmt.Figure()
fig.basemap(
    frame='a',
    projection=f"X10",
    region=[-10, 0, -10, 0]
)
fig.grdimage(da)
fig.show()

get us the following

so maybe we should allow it and only raise an error when coord_inc ==0?

[maxrjones]

so maybe we should allow it and only raise an error when coord_inc ==0?

Sounds good.

[seisman]

I'm reopening the issue because a similar fix should be done in the upstream GMT code.

@PaulWessel @meghanrjones I think the problem is the API function GMT_Create_Data doesn't check if the parameter wesn is reasonable (in this case, wesn=[0, 0, 0, 0] is passed).

[maxrjones]

I'm reopening the issue because a similar fix should be done in the upstream GMT code.

@PaulWessel @meghanrjones I think the problem is the API function GMT_Create_Data doesn't check if the parameter wesn is reasonable (in this case, wesn=[0, 0, 0, 0] is passed).

Thanks for the information. @PaulWessel, I would like to work on fixing this issue. I'll submit a PR in the GMT repository for your feedback and ping if help is needed. We should remember to run gmt.jl tests before merging anything.