grdimage shading fails with external matrix

[seisman]

The issue was first reported by @liamtoney in GenericMappingTools/pygmt#364. Using the PyGMT script in that issue, he got the following error:

grdimage [ERROR]: Passing zmax <= zmin prevents automatic CPT generation!
grdimage [ERROR]: Failed to read CPT geo.

After some debugging, it seems that GMT can correctly read the matrix (i.e., the earth relief data) and can plot it correctly if shading=False (i.e., no -I+d in the command line). However, if shading=True is used (i.e., grdimage -I+d), GMT reports the error above, possibly because the matrix in memory is modified/freed/re-initialized when calling grdgradient.

---

The issue can be partially reproduced using the following C codes. The C code below is modified from GMT's own test code testapi_matrix_plot.c. It reads a 2d matrix and plots the matrix using grdimage.

#include "gmt.h"

int main () {
	void *API = NULL;                /* The API control structure */
	struct GMT_MATRIX *M = NULL;    /* Structure to hold input matrix */
	char input[GMT_VF_LEN] = {""};	/* String to hold virtual input filename */
	char args[128] = {""};         	/* String to hold module command arguments */
	/* Initialize the GMT session */
	API = GMT_Create_Session ("test", 2U, GMT_SESSION_EXTERNAL, NULL);
	M = GMT_Read_Data (API, GMT_IS_MATRIX, GMT_IS_FILE, GMT_IS_SURFACE, GMT_READ_NORMAL, NULL, "2d-matrix.txt", NULL);
	/* Associate our matrix with a virtual file */
	GMT_Open_VirtualFile (API, GMT_IS_GRID|GMT_VIA_MATRIX, GMT_IS_SURFACE, GMT_IN, M, input);

	/* Prepare the module arguments */
	sprintf (args, "%s -JX6i -P -Baf -Cgeo", input);
	/* Call the grdimage module */
	GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, args);

	/* Close the virtual file */
	GMT_Close_VirtualFile (API, input);
	/* Destroy session */
	if (GMT_Destroy_Session (API)) return EXIT_FAILURE;
};

The data 2d-matrix.txt is the 60m earth relief data, converted from netCDF to matrix format using following Python code. The data is attached 2d-matrix.txt.

import numpy as np
import pygmt
er = pygmt.datasets.load_earth_relief("60m")
np.savetxt("2d-matrix.txt", np.flip(er.data, axis=0), fmt='%.1f')

The C code is compiled using

gcc testapi_matrix_plot.c $(gmt-config --cflags --libs)

The output image looks correct. The coordinates are incorrect, simply because I didn't pass the geographic information.

Changing the grdimage arguments to:

sprintf (args, "%s -JX6i -P -Baf -Cgeo -I+d", input);

to enable shading, then the output image is:

This is what I expect:

[joa-quim]

Hmmm, unless I'm missing something I think it's the user responsibility to pass arrays with order that GMT is expecting. I do that all the time in Julia and never had this problem. Well not counting the countless days working in a way that the memory layout can be set, and changed, to work with this.
Even netCDF can be TOPdown or BOTTOMup.

[seisman]

I didn't mean that the first figure is wrong. The flipping Y axis is simply because I didn't spend any time to re-sort the matrix. The first figure shows that the GMT_Read_Data can correctly read the matrix.

The main issue here is the second figure. After adding -I+d, the output figure is not what we expect.

[joa-quim]

Does it give the same thing if you create the ascii version with grd2xyz?

[seisman]

Sorry, I forgot to attach the data in my first comment 2d-matrix.txt.

[seisman]

Does it give the same thing if you create the ascii version with grd2xyz?

grd2xyz outputs XYZ triples or Z values, but GMT_Read_Data expects a matrix here (GMT_IS_MATRIX is passed to GMT_Read_Data).

[joa-quim]

Sorry but still fail to see that the problem is in GMT. As I said either in Mirone or Julia I pass grids back and forth and when such cases arise it was always my fault.
If the array has the correct the correct order, the illumination must have the same order. And speaking on orders, the header struct has a mem_layout member to deal with issues of this type.

[seisman]

@joa-quim It's not about the orders. I've updated the issue description. Hopefully it's more clear to read now.

[PaulWessel]

grdgradient does recycle the input grid and places gradients there except when we are not allowed:

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

So since In.file is a @GMTAPI@ memory reference it allocates a separate grid. So not immediately obvious why/where it overwrites the input grid. But clearly it does. @seisman, if you are able to walk through this in Xcode or similar you may be able to verify that it is not allocating a new grid. E.g., new_grid should be true in this case.

I can help with debug later but busy day with zoom with all Chairs and Provost. Shit-storm is brewing with Governor planning to take all vacant positions away (we have 5 with active searches all frozen...).

[seisman]

new_grid is true, but this line returns Grid_orig[0] with data=0.

gmt/src/grdimage.c

Line 891 in ddd829b

if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file[k], Grid_orig[k]) == NULL) { /* Get grid data */

[seisman]

The original Python script now crashes. It worked when I tested PR #3510, but I don't why now it crashes after merging #3510 to master.

import pygmt

fig = pygmt.Figure()

dem = pygmt.datasets.load_earth_relief('10m')

fig.grdimage(dem, region=[-180, 180, -90, 90], frame='af',
             projection='Cyl_stere/6i', cmap='geo', shading=True)

fig.show(method='external')

Here are the verbose messages reported by the above Python script.

log.txt

[PaulWessel]

Hm, I do not see any errors in that log.

[seisman]

It crashes at the end of the log.

[seisman]

Different from the C example code above, PyGMT calls following C API functions:

1. GMT_Create_Data
2. GMT_Put_Matrix
3. GMT_Open_VirtualFile

[seisman]

Re-open the issue because changes between GMT 6.1.0 and 6.1.1 breaks matrix shading again.

Here is a Python script to reproduce the issue (adapted from GenericMappingTools/pygmt#364):

import pygmt

fig = pygmt.Figure()

# shaded relief using @earth_relief_30m
fig.grdimage('@earth_relief_30m', region='g', frame=True, projection='Cyl_stere/6i', cmap='geo', shading=True)

fig.shift_origin(xshift='7i')

# shaded relief using xarray: DOES NOT WORK.
grid = pygmt.datasets.load_earth_relief('30m')
fig.grdimage(grid, region='g', frame=True, projection='Cyl_stere/6i', cmap='geo', shading=True)

fig.savefig("shading.png")

Actual output:

You can see that the right-side map doesn't have shading, although we use shading=True.

[seisman]

Ping @PaulWessel to the regression introduced in GMT 6.1.1, which breaks PyGMT's grdimage shading again.

[PaulWessel]

OK, I will fire up my PyGMT and Xcode tomorrow and have a look. I should remember just enough to at least tell where it goes wrong.

[PaulWessel]

OK, may need to refine my debug instructions for pygmt to be fool-proof: Got this far:

>>> import pygmt
>>> fig = pygmt.Figure()
>>> grid = pygmt.datasets.load_earth_relief('30m')
>>> fig.grdimage(grid, region='g', frame=True, projection='Cyl_stere/6i', cmap='geo', shading=True)
grdimage [ERROR]: File /Users/pwessel/opt/miniconda3/share/geo.cpt was not found
grdimage [ERROR]: Cannot find file /Users/pwessel/opt/miniconda3/share/geo.cpt
grdimage [ERROR]: File /Users/pwessel/opt/miniconda3/share/geo.cpt not found
[Session pygmt-session (5)]: Error returned from GMT API: GMT_FILE_NOT_FOUND (16)
grdimage [ERROR]: Failed to read CPT geo.

Do I need to set anything other than GMT_LIBRARY_PATH?

[seisman]

You're running your GMT debug build, right? For debug build, I believe you also need to define GMT_SHAREDIR.

[PaulWessel]

Thanks that works. It is a bit odd since the Xcode debug build does not install a share dir so I am pointing it to my regular non-xcode-debug build. So I have

GMT_LIBRARY_PATH=/Users/pwessel/GMTdev/gmt-dev/xbuild/src/Debug
GMT_SHAREDIR=/Users/pwessel/GMTdev/gmt-dev/dbuild/gmt6/share

Not a problem, but not sure to have the xcodebuild "install" the stuff under the xbuild tree.

[PaulWessel]

OK, so on first going-through everything looks fine: Grid is received, values look fine, they are passed to grdgradient which creates a new grid with intensities that are returned back, with values that look good. Then both grids are projected and that looks fine too. THen we build the rgb and I cannot see anything obviously wrong. I only stepped through a few nodes and the only thing I need to go back and do is to look at more values since the few I looked at had very small intensity values (0.0003 etc) which of course would not change the color. The projected intensity grid header has reasonable ranges though.
Is it easy for you to diff the two PNGs to see if they differ in subtle ways, as expected if the intensities are off by some unknown factor and close to 0?

[seisman]

Here is the diff image between the left and right plots:

For the right panel (xarray grid), I sometimes get a plot like this:

[PaulWessel]

Hm, so the diff image does seem to correlate a bit with data gradients, especially the Arctic continental shelf. Yet still a bit odd. And your white washing (I assume it is the same CPT?) would indicate a good amount of constant ambient light (positive intensities). OK, will contrast with the first command to see if values are radically different in the intensity grid.

[PaulWessel]

Not there yet, but (re-)learning that with the xarray there are no pad, of course, so the grdgradient call makes a new grid for the intensities with a pad of 1. This is returned to grdimage. The grid has no pad. Presumably something funny happens in this context. The left (file case) has pad = 2 throughout and life is good.

[PaulWessel]

More: Well ,Xcode is acting strange today. I step into a function, then when I step back out I only see assembly, not C. So that is a bit annoying. Anyhow, the new scheme with passing a common H header structure fails since the grid and the intensity grid has different paddings. So I certainly will need to pass HG and HI or something. So will do that first.

[PaulWessel]

@seisman, if I instead wanted to read in another (very small) grid entirely from netcdf into memory (like grid), is there support for that?

[seisman]

@seisman, if I instead wanted to read in another (very small) grid entirely from netcdf into memory (like grid), is there support for that?

Yes, you can use xarray.open_datatset to read a grid into memory. Here is an example.

I first generate a small grid using grdcut:
gmt grdcut @earth_relief_01d -R0/10/0/10 -Gsmall_grid.nc

then plot the grid in two different ways:

import pygmt
import xarray as xr

fig = pygmt.Figure()

fig.grdimage('small_grid.nc', frame=True, projection='Cyl_stere/6i', cmap='geo', shading=True)
fig.shift_origin(xshift='7i')

# shaded relief using xarray: DOES NOT WORK.
grid = xr.open_dataarray("small_grid.nc")
fig.grdimage(grid, frame=True, projection='Cyl_stere/6i', cmap='geo', shading=True)

fig.savefig("shading.png")