weiji14
commented
5 years ago Hmm, I've managed to reproduce this on my computer. Might be something to do with the virtualfile_from_grid mechanism, but I can't tell for sure. Will need to investigate.
Closed liamtoney closed 4 years ago
weiji14
commented
5 years ago Hmm, I've managed to reproduce this on my computer. Might be something to do with the virtualfile_from_grid mechanism, but I can't tell for sure. Will need to investigate.
liamtoney
commented
4 years ago My current work-around is to generate an illumination grid via grdgradient and apply that to my grdimage call, but it's pretty clunky.
liamtoney
commented
4 years ago 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')
ibrewster
commented
4 years ago @liamtoney For what it's worth, when I try your nested context manager"solution", I just get an error:
grdimage [ERROR]: Dimensions of intensity grid do not match that of the data grid!
Using the @ name, however, seems to work well - and runs MUCH faster, for some reason.
liamtoney
commented
4 years ago @liamtoney For what it's worth, when I try your nested context manager"solution", I just get an error:
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:
ibrewster
commented
4 years ago 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 shading=True not working (though, in my case, it actually crashes Python. Maybe Python version difference? Hmmm...)
liamtoney
commented
4 years ago @ibrewster I think it's an issue in regular GMT, see https://github.com/GenericMappingTools/gmt/issues/3408
seisman
commented
4 years ago It's confirmed to be a bug in GMT. Upstream PRs https://github.com/GenericMappingTools/gmt/pull/3510 and https://github.com/GenericMappingTools/gmt/pull/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.
liamtoney
commented
4 years ago I noticed another issue with shading, this time in grdview(), when making #502 — posting here because I think it's the same problem and should be fixed by the GMT upstream fix — @seisman could you confirm?
Here I construct the same grid using both grdmath (took me a bit!) and Python:
gmt 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.grdimport 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 INPUT = "data.grd", produces the correct result (at left) while INPUT = data produces an unexpected result (at right).
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()INPUT = "data.grd" |
INPUT = data |
|---|---|
 |
 |
seisman
commented
4 years ago @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.
liamtoney
commented
4 years ago I believe it's already fixed in upstream GMT.
Awesome, thanks for testing. I can edit the script in #502 to add I="+" once we support newer GMT.
liamtoney
commented
4 years ago Confirmed that both of the problematic examples above are fixed using PyGMT 1.2.0 and (the actual fix) GMT 6.1.0!
weiji14
commented
4 years ago We should add a test for this after GMT 6.1.0 is released.
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.
MarkWieczorek
commented
4 years ago I think that the docstring of grdimage needs to be updated. It is not clear what shading=True means (there is no similar option in gmt's grdimage). The default should be to provide an xarray grid, but this is not currently supported.
liamtoney
commented
4 years ago If the shading kwarg is not provided, then no shading is performed. If it is provided as shading=True, I think it should default to +d — see the GMT info here: https://docs.generic-mapping-tools.org/latest/grdimage.html#i. Currently, shading=True is the equivalent of -I with no args. According to the linked grdimage docs the default for a bare -I is no illumination. However, trying a simple GMT grdimage one-liner results in shading even if the bare -I is used. I.e., -I and -I+d return the same shaded image. (Side note: The gmt grdimage docs for the -I flag could be more clear on this.)
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 shading=True is equivalent to -I which (apparently) is equivalent to -I+d.
The default should be to provide an xarray grid, but this is not currently supported.
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.
MarkWieczorek
commented
4 years ago For reference, here is what the gmt man page says about the -I parameter (the web documentation is worded slightly differently):
usage: gmt grdimage ... -I[<intensgrid>|<value>|<modifiers]
and
-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). To derive from another grid than <grd_z>, give a filename as well.
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 +d uses the default values of option 3.
When the -I was first introduced, the only option was to specify a gridfile, and this usage persisted up until about gmt v5 when options 2 and 3 will introduced. This is why I just assumed that option 1 should be the default (i.e., I didn't know about the other 2 options until today!). However, +d might be a good fallback behaviour if -I is called with no other parameters.
weiji14
commented
4 years ago Looks like this got fixed in GMT upstream at https://github.com/GenericMappingTools/gmt/pull/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 https://github.com/conda-forge/gmt-feedstock/pull/108 so people can run conda install -c conda-forge/label/dev gmt to install a GMT development version to test things out now too :grin:
seisman
commented
4 years ago Closed by GMT upstream PR https://github.com/GenericMappingTools/gmt/pull/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)
``` grdimage [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter. grdimage [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt grdimage [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/cache grdimage [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server grdimage [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server/srtm1 grdimage [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server/srtm3 grdimage [DEBUG]: Got regular w/e/s/n for region (-180/180/-90/90) grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/cache grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server/srtm1 grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server/srtm3 grdimage [DEBUG]: Look for file @GMTAPI@-000000 in /Users/liamtoney/.gmt grdimage [DEBUG]: Look for file @GMTAPI@-000000 in /Users/liamtoney/.gmt/cache grdimage [DEBUG]: Look for file @GMTAPI@-000000 in /Users/liamtoney/.gmt/server grdimage [DEBUG]: Look for file @GMTAPI@-000000 in /Users/liamtoney/.gmt/server/srtm1 grdimage [DEBUG]: Look for file @GMTAPI@-000000 in /Users/liamtoney/.gmt/server/srtm3 grdimage [INFORMATION]: Read header from file @GMTAPI@-000000 grdimage [DEBUG]: Local file /Users/liamtoney/.gmt/server/gmt_hash_server.txt found grdimage [DEBUG]: File /Users/liamtoney/.gmt/server/gmt_hash_server.txt less than 24 hours old, refresh is premature. grdimage [DEBUG]: api_begin_io: Input resource access is now enabled [container] grdimage [DEBUG]: api_import_grid: Passed ID = 0 and mode = 1 grdimage [DEBUG]: Grid/Image dimensions imply w/e/s/n = -180/180/-90/90, inc = 0.166667/0.166667, gridline registration, n_layers = 1 grdimage [DEBUG]: Geographic input grid, longitudes span exactly 360 grdimage [INFORMATION]: Given domain implies x_inc = 0.166667 grdimage [INFORMATION]: Given domain implies y_inc = 0.166667 grdimage [DEBUG]: Grid is considered to have a 360-degree longitude range. grdimage [DEBUG]: Chosen boundary condition for all edges: geographic grdimage [DEBUG]: Geographic input grid, longitudes span exactly 360 grdimage [DEBUG]: Object ID 1 : Registered Grid Memory Reference 7fc1ced29590 as an Input resource with geometry Surface [n_objects = 2] grdimage [DEBUG]: Successfully created a new Grid container grdimage [DEBUG]: GMT_End_IO: Input resource access is now disabled grdimage [WARNING]: Spherical approximation used! grdimage [WARNING]: Central meridian not given, default to 0 grdimage [DEBUG]: Projected values in meters: -2.00151e+07 2.00151e+07 -1.2742e+07 1.2742e+07 grdimage [DEBUG]: Auto-frame interval for axis 0 item 0: d = 60 f = 15 grdimage [INFORMATION]: Auto-frame interval for x-axis (item 0): a60f15 grdimage [DEBUG]: Auto-frame interval for axis 1 item 0: d = 60 f = 15 grdimage [INFORMATION]: Auto-frame interval for y-axis (item 0): a60f15 grdimage [INFORMATION]: Map scale is 2626.65 km per cm or 1:2.62665e+08. grdimage [DEBUG]: Projected grid is non-orthogonal, nonlinear, or dpi was changed grdimage [INFORMATION]: Derive intensity grid from data grid grdimage [DEBUG]: Object ID 2 : Registered Grid Memory Copy 7fc1cec5a210 as an Output resource with geometry Surface [n_objects = 3] grdimage [DEBUG]: Successfully created a new Grid container grdimage [INFORMATION]: Calling grdgradient with args -G@GMTAPI@-000002 -A-45.0 -Nt1 -R-180/180/-90/90 --GMT_HISTORY=false @GMTAPI@-000000 grdimage [DEBUG]: GMT now running in modern mode [Session ID = 9323] grdimage (gmtlib_free_tmp_arrays): tried to free unallocated memory grdgradient [DEBUG]: History: Process -R-180/180/-90/90 grdgradient [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt grdgradient [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/cache grdgradient [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server grdgradient [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server/srtm1 grdgradient [DEBUG]: Look for file -180/180/-90/90 in /Users/liamtoney/.gmt/server/srtm3 grdgradient [DEBUG]: Got regular w/e/s/n for region (-180/180/-90/90) grdgradient [DEBUG]: Look for file -45.0 in /Users/liamtoney/.gmt grdgradient [DEBUG]: Look for file -45.0 in /Users/liamtoney/.gmt/cache grdgradient [DEBUG]: Look for file -45.0 in /Users/liamtoney/.gmt/server grdgradient [DEBUG]: Look for file -45.0 in /Users/liamtoney/.gmt/server/srtm1 grdgradient [DEBUG]: Look for file -45.0 in /Users/liamtoney/.gmt/server/srtm3 grdgradient [DEBUG]: amplitude_set = Y offset_set = N sigma_set = N grdgradient [INFORMATION]: Processing input grid grdgradient [DEBUG]: api_begin_io: Input resource access is now enabled [container] grdgradient [DEBUG]: api_import_grid: Passed ID = 0 and mode = 1 grdgradient [DEBUG]: Grid/Image dimensions imply w/e/s/n = -180/180/-90/90, inc = 0.166667/0.166667, gridline registration, n_layers = 1 grdgradient [DEBUG]: Geographic input grid, longitudes span exactly 360 grdgradient [INFORMATION]: Given domain implies x_inc = 0.166667 grdgradient [INFORMATION]: Given domain implies y_inc = 0.166667 grdgradient [DEBUG]: Grid is considered to have a 360-degree longitude range. grdgradient [DEBUG]: Chosen boundary condition for all edges: geographic grdgradient [DEBUG]: Geographic input grid, longitudes span exactly 360 grdgradient [DEBUG]: Object ID 3 : Registered Grid Memory Reference 7fc1cc5eb9f0 as an Input resource with geometry Surface [n_objects = 4] grdgradient [DEBUG]: Successfully created a new Grid container grdgradient [DEBUG]: GMT_End_IO: Input resource access is now disabled grdgradient [DEBUG]: api_begin_io: Input resource access is now enabled [container] grdgradient [DEBUG]: api_import_grid: Passed ID = 0 and mode = 2 grdgradient [INFORMATION]: Importing grid data from user memory location grdgradient [DEBUG]: Grid is considered to have a 360-degree longitude range. grdgradient [DEBUG]: Chosen boundary condition for all edges: geographic grdgradient [INFORMATION]: Set boundary condition for top edge: geographic grdgradient [INFORMATION]: Set boundary condition for bottom edge: geographic grdgradient [DEBUG]: GMT_End_IO: Input resource access is now disabled grdgradient [DEBUG]: Object ID 4 : Registered Grid Memory Reference 7fc1cc5d3850 as an Input resource with geometry Surface [n_objects = 5] grdgradient [DEBUG]: Successfully duplicated a Grid grdgradient [INFORMATION]: Min Mean Max sigma intensities: grdgradient [INFORMATION]: -0.35723270421 0.000249671630837 0.283594327547 0.0132078509689 grdgradient [DEBUG]: GMT_Write_Data: Writing Grid to memory object 2 from object 4 which transfers ownership grdgradient [DEBUG]: api_begin_io: Output resource access is now enabled [container] grdgradient [DEBUG]: api_export_data: Messenger dummy output container for object 2 [item 2] freed and set resource=data=NULL grdgradient [DEBUG]: api_export_grid: Passed ID = 2 and mode = 0 grdgradient [INFORMATION]: Referencing grid data to GMT_GRID memory location grdgradient [DEBUG]: GMT_End_IO: Output resource access is now disabled grdgradient [DEBUG]: gmtapi_garbage_collection: Destroying object: C=0 A=0 ID=3 W=Input F=Grid M=Memory Copy S=Used P=7fc1cc5eb9f0 N=(null) grdgradient [DEBUG]: GMTAPI_Garbage_Collection freed 1 memory objects grdgradient [DEBUG]: gmtapi_unregister_io: Unregistering object no 3 [n_objects = 4] grdgradient [DEBUG]: gmtapi_unregister_io: Unregistering object no 4 [n_objects = 3] grdgradient (gmtlib_free_tmp_arrays): tried to free unallocated memory grdimage [DEBUG]: Running in PS mode modern grdimage [DEBUG]: Use PS filename /Users/liamtoney/.gmt/sessions/gmt6.9323/gmt_8.ps- grdimage [DEBUG]: Create hidden PS file /Users/liamtoney/.gmt/sessions/gmt6.9323/gmt_8.ps- grdimage [DEBUG]: Got session name as pygmt-session and default graphics formats as pdf grdimage [INFORMATION]: Allocate and read data from file @GMTAPI@-000000 grdimage [DEBUG]: api_begin_io: Input resource access is now enabled [container] grdimage [DEBUG]: api_import_grid: Passed ID = 1 and mode = 2 grdimage [INFORMATION]: Duplicating grid data from GMT_GRID memory location grdimage [DEBUG]: GMT_End_IO: Input resource access is now disabled grdimage [INFORMATION]: Allocates memory and read intensity file grdimage [INFORMATION]: project grid files grdimage [DEBUG]: Object ID 5 : Registered Grid Memory Reference 7fc1cc765be0 as an Input resource with geometry Surface [n_objects = 4] grdimage [DEBUG]: Successfully duplicated a Grid grdimage [DEBUG]: gmt_project_init: IN: Inc [0/0] n_columns/n_rows [2161/1081] dpi = 0 offset = 0 grdimage [DEBUG]: gmt_project_init: OUT: Inc [0/0] n_columns/n_rows [2161/1081] dpi = 0 offset = 0 grdimage [INFORMATION]: Grid projection from size 2161x1081 to 2161x1081 grdimage [DEBUG]: Successfully added data array to previously registered Grid container grdimage [DEBUG]: gmt_grd_project: In [-180/180/-90/90] and out [0/6/0/3.81971863421] grdimage [DEBUG]: GMT_Destroy_Data: freed memory for a Grid for object 1 grdimage [DEBUG]: gmtapi_unregister_io: Unregistering object no 1 [n_objects = 3] grdimage [DEBUG]: gmtapi_unregister_io: Object no 1 has non-NULL resource pointer grdimage [DEBUG]: Object ID 6 : Registered Grid Memory Reference 7fc1cc73a3e0 as an Input resource with geometry Surface [n_objects = 4] grdimage [DEBUG]: Successfully duplicated a Grid grdimage [DEBUG]: gmt_project_init: IN: Inc [0/0] n_columns/n_rows [2161/1081] dpi = 0 offset = 0 grdimage [DEBUG]: gmt_project_init: OUT: Inc [0/0] n_columns/n_rows [2161/1081] dpi = 0 offset = 0 grdimage [INFORMATION]: Grid projection from size 2161x1081 to 2161x1081 grdimage [DEBUG]: Successfully added data array to previously registered Grid container grdimage [DEBUG]: gmt_grd_project: In [-180/180/-90/90] and out [0/6/0/3.81971863421] grdimage [WARNING]: gmt_grd_project: Output grid extrema [-1.89669/1.87961] exceed extrema of input grid [-0.97649/0.970346] due to resampling grdimage [WARNING]: gmt_grd_project: Output grid clipped to input grid extrema grdimage [DEBUG]: GMT_Destroy_Data: freed memory for a Grid for object 2 grdimage [DEBUG]: gmtapi_unregister_io: Unregistering object no 2 [n_objects = 3] grdimage [DEBUG]: gmtapi_unregister_io: Object no 2 has non-NULL resource pointer grdimage [INFORMATION]: Given grid header, select default CPT to be turbo grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/cache grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server/srtm1 grdimage [DEBUG]: Look for file geo in /Users/liamtoney/.gmt/server/srtm3 grdimage [ERROR]: Passing max <= zmin prevents automatic CPT generation! grdimage [ERROR]: Failed to read CPT geo. grdimage [DEBUG]: gmtapi_garbage_collection: Destroying object: C=0 A=0 ID=5 W=Input F=Grid M=Memory Reference S=Unused P=7fc1cc765be0 N=(null) grdimage [DEBUG]: gmtapi_garbage_collection: Destroying object: C=0 A=0 ID=6 W=Input F=Grid M=Memory Reference S=Unused P=7fc1cc73a3e0 N=(null) grdimage [DEBUG]: GMTAPI_Garbage_Collection freed 2 memory objects grdimage [DEBUG]: gmtapi_unregister_io: Unregistering object no 5 [n_objects = 2] grdimage [DEBUG]: gmtapi_unregister_io: Unregistering object no 6 [n_objects = 1] grdimage (gmtlib_free_tmp_arrays): tried to free unallocated memory Error: /syntaxerror in --%ztokenexec_continue-- Operand stack: PSL_draw_path_lines Execution stack: %interp_exit .runexec2 --nostringval-- --nostringval-- --nostringval-- 2 %stopped_push --nostringval-- --nostringval-- --nostringval-- false 1 %stopped_push 1999 1 3 %oparray_pop 1998 1 3 %oparray_pop 1982 1 3 %oparray_pop 1868 1 3 %oparray_pop --nostringval-- %errorexec_pop .runexec2 --nostringval-- --nostringval-- --nostringval-- 2 %stopped_push --nostringval-- Dictionary stack: --dict:981/1684(ro)(G)-- --dict:0/20(G)-- --dict:78/200(L)-- --dict:136/250(L)-- Current allocation mode is local Last OS error: No such file or directory psconvert [ERROR]: System call [gs -q -dSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL_no_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true '/Users/liamtoney/.gmt/sessions/gmt6.9323/gmt_8.ps-' 2> '/Users/liamtoney/.gmt/sessions/gmt6.9323/psconvert_9323c.bb'] returned error 256. --------------------------------------------------------------------------- GMTCLibError Traceback (most recent call last)
System information
conda listbelow:output of conda list