Several problems in grdgradient

[PaulWessel]

grdgradient computes derivatives and requires at least a one row/col padding for the algorithm to work: We do the finite differencing on a cell defined by 4 nodes and store the result in the top left node which is not used in further calculations. At the end, we restore the padding so the output grid is correctly aligned. For files (automatically having a pad of 2) this works fine.

When grdgradient is called via GMT_Call_Module with a memory grid that has no padding (#3408) we detect that case and call

new_grid = gmt_set_outgrid (GMT, Ctrl->In.file, separate, 1, Surf, &Out); /* true if input is a read-only array */

This function returns true in this case and allocates a new grid (Out) that (given the arg = 1) will have the required minimum of one row/col padding. We then place computed gradients into Out->data[index] so that Surf (read-only) is not touched. But wait a minute! The loop we set up is using the Surf->header to get the index into the Surf array, but in this case Surf has no padding. Note we always do this:

	mx = Surf->header->mx;	/* Need a signed mx for p[3] in line below */
	p[0] = 1;	p[1] = -1;	p[2] = mx;	p[3] = -mx;

so p[1] and p[3] are always negative relative node changes. Then in the loop we do this

for (n = 0, bad = false; !bad && n < 4; n++) if (gmt_M_is_fnan (Surf->data[ij+p[n]])) bad = true;

Since Surf->data has no pad and ij starts at 0, we end up with negative array indices. Not sure why it is not crashing but it clearly points to nowhere. This explains part of the problem in #3408 but regardless this is a serious bug.

Another issues is -S. With -S we create another output grid Slope by duplicating Surf. We then use the same index as Out to assign output. But this will be wrong if Surf has no pad. Now I allocate Slope using the same padding as the output grid.

Finally, ensuring that the grid has at least one row/col is useless since gmt_set_outgrid may in fact need to call gmt_grd_BC_set which requires two rows/cols. So now we call it with two.

For clarity, this PR replaces "Surf" with "In", and "Out" with "Grid" and uses Grid in all those places since it is this array we get values from and then we always place the result in the top-left node of the cells. At the end we reinstate the pad. Doing this and running all the test gave no errors, and it also fixes #3408.

[seisman]

@PaulWessel This PR works well for most cases but crashes for a constant intensity for both CLI and PyGMT.

Here is the CLI crash:

$ gmt grdimage @earth_relief_01d -Baf -Cgeo -I0.5 -JCyl_stere/6i -R-180/180/-90/90 -pdf map
ERROR: Caught signal number 11 (Segmentation fault: 11) at
0   libgmt.6.dylib                      0x00000001089d5100 GMT_grdimage + 19264
1   ???                                 0x0000000000000000 0x0 + 0
Stack backtrace:
0   libgmt.6.dylib                      0x0000000108748332 sig_handler + 594
1   libsystem_platform.dylib            0x00007fff6cabd5fd _sigtramp + 29
2   ???                                 0x00007ffc4600ba00 0x0 + 140721482938880
3   libgmt.6.dylib                      0x0000000108769bbc GMT_Call_Module + 2156
4   gmt                                 0x000000010872f593 main + 1459
5   libdyld.dylib                       0x00007fff6c8c4cc9 start + 1
6   ???                                 0x000000000000000a 0x0 + 10

[PaulWessel]

All those

struct GMT_GRID_HEADER *H_i = Conf->Intens->header;

needs to be:

struct GMT_GRID_HEADER *H_i = (Conf->Intens) ? Conf->Intens->header : NULL;

Heading out if you want to do this.

[PaulWessel]

I'm back, unless @seisman started on this I can do it.

[seisman]

I'm not working on it.