Allow the use of a separate file for shading (grdimage -I)

[MarkWieczorek]

The gmt command grdimage has an option -I that adds shading to another image by modulating the intensity. Historically, (up until gmt v5), this command has only taken as an option another grd filename. One common use of this feature was to plot something like a gravity, the geoid, or magnetic field in color, and then apply shading using a topography file to provide geologic context.

Since gmt 5, there has been the option of using the input grd file along with grdgradient to automatically generate an intensity file that would be used with the -I command. I believe that pygmt only accepts the option shading=True which automatically generates such an intensity file using gmt defined default illumination angle parameters.

I am asking that pygmt fully implement the grdimage -I option. This is the only feature that I am currently missing in pygmt, and which is forcing me to use both gmt and pygmt separately to generate publication quality images for journal articles.

For reference, here is the relevant part of the gmt grdimage man page:

	-I Apply directional illumination. Append name of intensity grid file.
	   For a constant intensity (i.e., change the ambient light), append a value.
	   To derive intensities from <grd_z> instead, append +a<azim> [-45], +n<method> [t1], and +m<ambient> [0]
	   or use -I+d to accept the default values (see grdgradient for details).

[weiji14]

Hi @MarkWieczorek, I agree that this will be a useful feature. Do you have an example for this/a map illustration? Even just the pure GMT command will do (best if we can access the data too).

[MarkWieczorek]

Well, here is what it should look like (Earth geoid with topography illumination). Sorry for the rainbow colormap, this is a few years old:

I could look through my files, but the scripts will most likely be for gmt 5, and the grd files are pretty big. The only peculiarity with specifying a -I file is that the grid size and grid registration need to be exactly the same. If you really need something to test, I could upload them somewhere. But, the easiest test would be to just use the same grid for the colormap as for -I and see if you get the same thing as -I=True.

[weiji14]

I was thinking more along the lines of a GMT Remote Data-Set (i.e. the @ files), but not sure if there's one that would be suitable? Could you please post the code too for generating this map? Sorry for the fuss, just that I'm not familiar with shading.

[MarkWieczorek]

Here is something for a different image

grdimage -R-0/360/-90/90 thick.grd -I/Users/lunokhod/Mars/Topography/Mars_16_srm.grd -Csingle.cpt -JW-100/13 -Bg30WSen -Y10 -K -V > output.ps

In essence, if you are creating an image with grdimage, you just add the option -Ifile.grd

[weiji14]

Here is something for a different image

grdimage -R-0/360/-90/90 thick.grd -I/Users/lunokhod/Mars/Topography/Mars_16_srm.grd -Csingle.cpt -JW-100/13 -Bg30WSen -Y10 -K -V > output.ps

In essence, if you are creating an image with grdimage, you just add the option -Ifile.grd

Cool, thanks for the example. This sounds very similar to the drapegrid/-G option in pygmt.grdview, and it could perhaps be implemented similarly:

pygmt/pygmt/base_plotting.py

Lines 433 to 444 in d3e131a

	with contextlib.ExitStack() as stack:
	fname = stack.enter_context(file_context)
	if "G" in kwargs:
	drapegrid = kwargs["G"]
	if data_kind(drapegrid) in ("file", "grid"):
	if data_kind(drapegrid) == "grid":
	drape_context = lib.virtualfile_from_grid(drapegrid)
	drapefile = stack.enter_context(drape_context)
	kwargs["G"] = drapefile
	else:
	raise GMTInvalidInput(
	f"Unrecognized data type for drapegrid: {type(drapegrid)}"

I'll take a closer look over the weekend (getting late for me here), but you're more than welcome to open up a Pull Request to speed up the process 😄.

[seisman]

I believe that pygmt only accepts the option shading=True which automatically generates such an intensity file using gmt defined default illumination angle parameters.

I don't test it, but I believe you can already use shading="/Users/lunokhod/Mars/Topography/Mars_16_srm.grd" to speify the intensity grid. Please let us know if it works for you.

Of course, there are more things to do:

1. the shading is not well documented in the grdimage function
2. currently we can't pass an xarray grid to shading.

[MarkWieczorek]

I'll test (1) monday, but (2) is what I really want!

[seisman]

Yes, (2) is definitely a missing feature!

[weiji14]

I don't test it, but I believe you can already use shading="/Users/lunokhod/Mars/Topography/Mars_16_srm.grd" to speify the intensity grid. Please let us know if it works for you.

I can confirm that this works, here's a very minimal working example:

import pygmt

fig = pygmt.Figure()
fig.grdimage(grid="@earth_relief_01d_g", shading="@earth_relief_01d_g+d")
fig.savefig("shading_with_intensfile.png")
fig.show()

However, shading an xarray grid with an intensfile doesn't work, likely due to same issue at #364/#581:

fig = pygmt.Figure()
grid = pygmt.datasets.load_earth_relief(registration="gridline")
fig.grdimage(grid=grid, cmap="geo", shading="@earth_relief_01d_g+d")
fig.savefig("xarray_shading_with_intensfile.png")
fig.show()

Of course, there are more things to do:

1. the shading is not well documented in the grdimage function

2. currently we can't pass an xarray grid to shading.

Agree with both points, we can better document grdimage and allow for xarray grids to shading in the same PR. But note that for point 2 (shading with xarray grids), whatever we do won't work properly until we fix #364/#581 (an upstream GMT issue).

[seisman]

It seems that the upstream GMT PR GenericMappingTools/gmt#4328 also fixes this issue.

Running @weiji14's script above gives me a plot which looks correct (with the GMT master branch version: 6.2.0_8a4275f_2020.10.14):

But we need some more tests.

[MarkWieczorek]

I finally got around to testing this, and grdimage with shading works fine for me as long as I specify a pre-existing filename. For my pipeline, a hack solution is to save an xarray DataArray as a netcdf file, and then to use this filename with the shading option.

Obviously, it would be better if we could just specify the xarray DataArray directly in the shading option.

[seisman]

@MarkWieczorek I've implemented the requested feature in #750 and it seems working well (I only did a few simple tests and eye-checking the images). It would be good if you can test it and give some feedbacks.

[seisman]

The feature request has been implemented in #750 and now merged in master. Now grdimage can accept a netCDF grid or a xarray.DataArray as shading.

However, the grdimage shading with a xarray.DataArray has known upstream bugs (see #750 (comment), #750 (comment), and GenericMappingTools/gmt#5294).

I'm reopening this issue to better track the upstream bug.

[weiji14]

However, the grdimage shading with a xarray.DataArray has known upstream bugs (see #750 (comment), #750 (comment), and GenericMappingTools/gmt#5294).

I'm reopening this issue to better track the upstream bug.

The last bits of the upstream CPT stretching issue (GenericMappingTools/gmt#5294) should be fixed with GenericMappingTools/gmt#5947 and GenericMappingTools/gmt#5948 (thanks Meghan and Paul!) and will be available in GMT 6.3.

I've opened a GMT conda-forge dev build at conda-forge/gmt-feedstock#173 which will contain these patches, and allow us to test things easily using conda install -c conda-forge/label/dev gmt. So we should be able to close this issue once PyGMT bumps the minimum GMT version to 6.3, and the xfail is removed from the test below:

pygmt/pygmt/tests/test_grdimage.py

Lines 120 to 125 in 5dc135b

	@pytest.mark.xfail(
	reason="Incorrect scaling of geo CPT on xarray.DataArray grdimage plot."
	"See https://github.com/GenericMappingTools/gmt/issues/5294",
	)
	@check_figures_equal()
	def test_grdimage_grid_and_shading_with_xarray(grid, xrgrid):

The test is currently XPASS at https://github.com/GenericMappingTools/pygmt/runs/4111894254?check_suite_focus=true#step:15:295 😁