Allow passing xr.DataArray as shading to grdimage

[seisman]

Description of proposed changes

Following the comment #618 (comment) and the similar implementation in grdview function, this PR allows passing xr.DataArray shading to grdimage.

Here are two examples:

The shading.grd file is generated using GMT command, which is not available in PyGMT yet:

gmt grdgradient @earth_relief_01d_g -Gshading.grd -A-45 -Nt1+a0

importpygmtimportxarrayasxr# reading shading grid as xr.DataArrayshading=xr.open_dataarray("shading.grd")
# grid is from a file, shading is from xr.DataArrayfig=pygmt.Figure()
fig.grdimage(grid="@earth_relief_01d_g", shading=shading)
fig.show()
# both grid and shading are from xr.DataArrayfig=pygmt.Figure()
grid=pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, cmap="geo", shading=shading)
fig.show()

Fixes#618.

ISSUES

If shading is given as a grid file, GMT can automatically derive the gradient/intensity/shading by adding +d to the grid file name, e.g., fig.grdimage(grid="@earth_relief_01d_g", shading="@earth_relief_01d_g+d"). However, it's impossible when shading is given as a xr.DataArray. Perhaps we could add a new argument (e.g., shading_parameters="+d") and automatically append it to the shading file name or virtual file name.

TODO

• Add tests

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to docstrings or tutorials.

Notes

• You can write /format in the first line of a comment to lint the code automatically

[MarkWieczorek]

I will test this later in the weekend. As for the "+d" option, here is what I have been doing in pyshtools. I have two optional parameters: shading and shading_azimuth. (There is not enough info in the grdimage -I documentation for me to figure out what the +nt1+m0 means, so I am just keeping this as a default. I am not sure if these two parameters are of interest for shading.)

shading=True uses the base grid to create the shading, shading=filename sends the grd filename to GMT as is, and shading=xarray creates a temporary grd file that gets sent to GMT as is. This PR should avoid writing the temporary file to disk.

Using the input value of shading_azimuth, I create the string "+a{:}+nt1+m0".format(shading_azimuth), which gets appended to whatever shading is. If shading_azimuth is None, then nothing gets appended to shading, and the file is assumed to already be a gradient file.

As you mention, as opposed to using the "+d" option, we could just instead create the gradient map in pygmt using grdgradient (when it is wrapped). I'm not sure if this would be preferable to GMT creating the gradient or not.

[seisman]

There is not enough info in the grdimage -I documentation for me to figure out what the +nt1+m0 means

Yes, the grdimage -I documentation could be improved. +nt1+m0 is equivalent to -Nt1+a0 in grdgradient (https://docs.generic-mapping-tools.org/dev/grdgradient.html#n).

shading=xarray creates a temporary grd file that gets sent to GMT as is. This PR should avoid writing the temporary file to disk.

Yes, this PR uses the xarary in memory directly, so no temporary file on disk.

[weiji14]

The code implementation is technically correct already in pygmt/src/grdimage.py, just that the xarray.DataArray grdimage plot appears different to that from a NetCDF grid. What do you think about merging this PR into master (after adding an xfail to the test) so that can Paul test things easier in GenericMappingTools/gmt#5294?

[weiji14]

Btw, I tried using makecpt to get both plots to use the same colormap scale, but the images still differ, appears to be some slight offset in the left-right direction. So it's not just a matter of CPT scaling as reported in GenericMappingTools/gmt#5294, but maybe 2 problems.

Suggested change

	fig_ref.grdimage(
	grid="@earth_relief_01d_g", region="GL", cmap="geo", shading=xrgrid, verbose="i"
	)
	fig_ref.colorbar()
	fig_test.grdimage(grid=grid, region="GL", cmap="geo", shading=xrgrid, verbose="i")
	fig_test.colorbar()
	makecpt(cmap="geo", series=(-3900, 3200))
	fig_ref.grdimage(
	grid="@earth_relief_01d_g", region="GL", cmap=True, shading=xrgrid, verbose="d"
	)
	fig_ref.colorbar()
	fig_test.grdimage(grid=grid, region="GL", cmap=True, shading=xrgrid, verbose="d")
	fig_test.colorbar()

Expected (correct)	Generated (wrong)	Diff

[seisman]

The code implementation is technically correct already in pygmt/src/grdimage.py, just that the xarray.DataArray grdimage plot appears different to that from a NetCDF grid. What do you think about merging this PR into master (after adding an xfail to the test) so that can Paul test things easier in GenericMappingTools/gmt#5294?

Sounds good to me.

[weiji14]

Good work! I think we can merge after the tests pass so that the debugging work at GenericMappingTools/gmt#5294 can proceed.