grdimage shading fails for xarray input #364
Comments
|
Hmm, I've managed to reproduce this on my computer. Might be something to do with the |
|
My current work-around is to generate an illumination grid via |
|
Here it is, in all of its twice-nested-context-manager glory: import pygmt
fig = pygmt.Figure()
dem = pygmt.datasets.load_earth_relief('10m')
with pygmt.helpers.GMTTempFile() as tmp_grd:
with pygmt.clib.Session() as session:
with session.virtualfile_from_grid(dem) as dem_file:
session.call_module('grdgradient', f'{dem_file} -A-45 -Nt1- -G{tmp_grd.name}')
fig.grdimage(dem, region=[-180, 180, -90, 90], frame='af',
projection='Cyl_stere/6i', cmap='geo', shading=tmp_grd.name)
fig.show(method='external') |
|
@liamtoney For what it's worth, when I try your nested context manager"solution", I just get an error:
Using the |
Are you running the exact same code as above? Perhaps using a different grid make yours not work? I just tested the above code again on the GMT Binder, and it worked fine. You can try it here: |
|
Smaller grid (tried at 01m and 30s), smaller region, and different projection, but otherwise copied-and-pasted your code, yes. Of course, this is still just a workaround for |
|
@ibrewster I think it's an issue in regular GMT, see GenericMappingTools/gmt#3408 |
|
It's confirmed to be a bug in GMT. Upstream PRs GenericMappingTools/gmt#3510 and GenericMappingTools/gmt#3512 fixed the bug and were already merged into GMT's master branch. If you build the latest GMT master branch, running @liamtoney's example script in the top post will give you a beautiful earth relief map with shadings.
For users who don't want to manually build the GMT source codes, you can wait for the upcoming GMT 6.1.0, which is planned to be released before July 1st. NOTE for PyGMT maintainers: We should add a test for this after GMT 6.1.0 is released. |
|
I noticed another issue with shading, this time in Here I construct the same grid using both Pure GMTgmt grdmath -R-5/5/-5/5 -I0.05 X Y HYPOT 0.5 0.5 POW MUL -0.2 MUL EXP -20 MUL 2 PI MUL X MUL COS 2 PI MUL Y MUL COS ADD 0.5 MUL EXP SUB E ADD 20 ADD = data.grdPythonimport numpy as np
import xarray
def ackley(x, y):
return (
-20 * np.exp(-0.2 * np.sqrt(0.5 * (x ** 2 + y ** 2)))
- np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y)))
+ np.exp(1)
+ 20
)
INC = 0.05
x = np.arange(-5, 5 + INC, INC)
y = np.arange(-5, 5 + INC, INC)
data = xarray.DataArray(ackley(*np.meshgrid(x, y)), coords=(x, y))Then the following script, with import pygmt
fig = pygmt.Figure()
SCALE = 0.2 # [inches]
fig.grdview(
INPUT,
frame=["a5f1", "za5f1"],
projection=f"x{SCALE}i",
zscale=f"{SCALE}i",
surftype="s",
cmap="roma",
perspective="135/30",
I="+",
)
fig.show()
|
|
@liamtoney It's a nice figure and a nice test! I can reproduce your issue with GMT 6.0.0, but CAN'T reproduce it with GMT master. I believe it's already fixed in upstream GMT. |
Awesome, thanks for testing. I can edit the script in #502 to add |
|
Confirmed that both of the problematic examples above are fixed using PyGMT 1.2.0 and (the actual fix) GMT 6.1.0! |
Ok, GMT 6.1.0 is released and we've bumped the master branch of PyGMT to use it (thanks @seisman!). If someone wants to learn a bit about writing unit tests, this might an easy one to have a go at. |
|
I think that the docstring of grdimage needs to be updated. It is not clear what |
|
If the In any case, I agree that the docstring for grdimage needs to be updated. The docstring for grdview (dev docs) is here: https://www.pygmt.org/dev/api/generated/pygmt.Figure.grdview.html#pygmt.Figure.grdview. We could use the same for grdimage and grdview, but add in (for both) that
I don't understand what you mean by default in this context, but I do think this should be supported. Probably good to make a new feature request for this though. |
|
For reference, here is what the gmt man page says about the
and
Thus, option 1 is to use a gridfile, option 2 is for a constant value, and option 3 is to derive intensities using grdgradient. If you chose option 3, then When the |
|
Looks like this got fixed in GMT upstream at GenericMappingTools/gmt#4328, I see the xarray shading test we added in #581 is now an xpass instead of xfail. So anyone who builds GMT from source (or wait for GMT 6.2.0 in Feb/Mar 2021?) should not face this bug anymore. In the meantime, I've made a PR at conda-forge/gmt-feedstock#108 so people can run |
The xarray shading issue in #364 was fixed by upstream GMT in GenericMappingTools/gmt#4328. This PR updates the pytest xfail condition so that the xarray shading test is expected to xfail only for GMT<=6.1.1.
|
Closed by GMT upstream PR GenericMappingTools/gmt#4328. |
Description of the problem
Passing
shading=Truein a simplefig.grdimage()call produces an error if the input grid is anxarray.DataArray. The error does not appear when I use the special "@" filename.Full code that generated the error
This works, though:
Error message summary
Full error message (debug verbosity)
System information
conda listbelow:output of conda list
The text was updated successfully, but these errors were encountered: