Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow passing xr.DataArray as shading to grdimage #750

Merged
merged 23 commits into from
Jun 10, 2021

Conversation

seisman
Copy link
Member

@seisman seisman commented Dec 19, 2020

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
import pygmt
import xarray as xr

# reading shading grid as xr.DataArray
shading = xr.open_dataarray("shading.grd")

# grid is from a file, shading is from xr.DataArray
fig = pygmt.Figure()
fig.grdimage(grid="@earth_relief_01d_g", shading=shading)
fig.show()

# both grid and shading are from xr.DataArray
fig = 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

@seisman seisman added enhancement Improving an existing feature feature Brand new feature and removed enhancement Improving an existing feature labels Dec 19, 2020
@seisman seisman changed the title Pass xr.DataArray shading to grdimage Allow passing xr.DataArray as shading to grdimage Dec 19, 2020
@seisman seisman added this to the 0.3.0 milestone Dec 19, 2020
seisman added a commit that referenced this pull request Dec 19, 2020
When I worked on PR #750 based on the implemention of the `grdview` function,
I found that the `grdview` codes are a little difficult to understand.
The main reason is that `contextlib.ExitStack()` is new to me.
Another reason is that these two codes `fname = stack.enter_context(file_context)`
and `arg_str = " ".join([fname, build_arg_string(kwargs)])` are
separated by the long `if` block.
@MarkWieczorek
Copy link
Contributor

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 added a commit that referenced this pull request Dec 19, 2020
When I worked on PR #750 based on the implemention of the `grdview` function,
I found that the `grdview` codes are a little difficult to understand.
The main reason is that `contextlib.ExitStack()` is new to me.
Another reason is that these two codes `fname = stack.enter_context(file_context)`
and `arg_str = " ".join([fname, build_arg_string(kwargs)])` are
separated by the long `if` block.
@seisman
Copy link
Member Author

seisman commented Dec 19, 2020

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.

Copy link
Member

@weiji14 weiji14 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

pygmt/tests/test_grdimage.py Show resolved Hide resolved
Comment on lines +127 to +132
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()
Copy link
Member

@weiji14 weiji14 Jun 9, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
test_grdimage_grid_and_shading_with_xarray png -expected test_grdimage_grid_and_shading_with_xarray png test_grdimage_grid_and_shading_with_xarray png -failed-diff

@seisman
Copy link
Member Author

seisman commented Jun 9, 2021

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.

pygmt/src/grdimage.py Outdated Show resolved Hide resolved
pygmt/src/grdimage.py Outdated Show resolved Hide resolved
Co-authored-by: Wei Ji <[email protected]>
@seisman seisman marked this pull request as ready for review June 10, 2021 02:06
Copy link
Member

@weiji14 weiji14 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@seisman seisman merged commit d26fe0a into master Jun 10, 2021
@seisman seisman deleted the grdimage-xarray-shading branch June 10, 2021 02:34
sixy6e pushed a commit to sixy6e/pygmt that referenced this pull request Dec 21, 2022
…s#750)

* Pass xr.DataArray shading to grdimage
* Test grdimage with xarray.DataArray input to both grid and shading
* Update grdimage shading docstring to mention xarray.DataArray input
* Mark the grdimage shading test as xfail [skip ci]

Co-authored-by: Wei Ji <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Improving an existing feature
Projects
No open projects
Release v0.3.x
  
Awaiting triage
Development

Successfully merging this pull request may close these issues.

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