Error pygmt GMTCLibError: Module 'grdimage' failed with status code 78

Saren · August 30, 2020, 10:50am

Hello, I have a problem while creating a map via grdimage and I don’t know how to fix it…

GMTCLibError: Module 'grdimage' failed with status code 78:
grdblend [ERROR]: Failed to remove C:\Users\Maureen\AppData\Local\Temp/grdblend_resampled_13784_2.nc! [remove error: Permission denied]
grdblend [ERROR]: Failed to delete file C:\Users\Maureen\AppData\Local\Temp/grdblend_resampled_13784_2.nc
grdimage [ERROR]: Dimensions of intensity grid do not match that of the data grid!

This is the code I wrote:

import pygmt

region=[3.970,4.270,44.000,44.250]
fig=pygmt.Figure()
fig.basemap(region=region, projection="M8i", frame=True)
fig.grdimage("@earth_relief_03s", cmap='geo', shading=True)
fig.coast(rivers='a/0.5,skyblue', water='skyblue')
fig.coast(borders='1/1,black')
fig.savefig('ales.png')

weiji14 · August 30, 2020, 10:09pm

Hi @Saren,

This is most likely due to the ‘shading’ setting, and it looks similar to the issue reported at https://github.com/GenericMappingTools/pygmt/issues/364 which should have been fixed with GMT 6.1.0, and it might work with PyGMT v0.1.2. However, I can reproduce this error on https://github.com/GenericMappingTools/try-gmt which uses these latest versions. I’ll do a bit more investigation to see where the problem is.

In the meantime, could you let us know what version of PyGMT/GMT you are using? Try doing pygmt.show_versions(), or if that doesn’t work, use pygmt.__version__() inside Python and gmt --version on the command-line.

weiji14 · August 30, 2020, 10:39pm

Ok, I think this might be an upstream GMT bug (ping @admins) and it doesn’t look to be fixed in the latest GMT 6.1 branch. This is the equivalent GMT bash script:

#!/usr/bin/env bash
gmt begin ales png
    gmt basemap -R3.970/4.270/44.000/44.250 -JM8i -B
    gmt grdimage @earth_relief_03s -Cgeo -I -V
    gmt coast -Ia/0.5,skyblue -Sskyblue
    gmt coast -N1/1,black
gmt end show

and the error message:

gmt [INFORMATION]: Reading grid from file srtm_tiles.nc
gmt [INFORMATION]: Set boundary condition for top    edge: natural
gmt [INFORMATION]: Set boundary condition for bottom edge: natural
grdblend [INFORMATION]: Reading Data Table from file ~/.gmt/sessions/gmt_session.27776/=tiled_80_PO.000000
grdblend [INFORMATION]: Given domain implies x_inc = 0.000833333
grdblend [INFORMATION]: Given domain implies y_inc = 0.000833333
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for latitude from 0.00416666666666666834 to 0.00416666666666666661
grdblend [INFORMATION]: Blend file ~/.gmt/server/earth/earth_relief/earth_relief_03s_g/N44E003.earth_relief_03s_g.nc in 3/4/44/45 with normal weight 1 [-900-300]
grdblend [INFORMATION]: Blend file ~/.gmt/server/earth/earth_relief/earth_relief_03s_g/N44E004.earth_relief_03s_g.nc in 4/5/44/45 with normal weight 1 [-900-300]
grdblend [WARNING]: File @N40E000.earth_relief_15s_p.nc has different increments (0.00416666666667/0.00416666666667) than the output grid (0.000833333333333/0.000833333333333) - must resample
grdblend [WARNING]: File @N40E000.earth_relief_15s_p.nc coordinates are phase-shifted w.r.t. the output grid - must resample
grdblend [INFORMATION]: Resample @N40E000.earth_relief_15s_p.nc via grdsample @N40E000.earth_relief_15s_p.nc -rg -I0.000833333333333/0.000833333333333 -R3.96666666667/4.27083333333/44/44.25 -G/tmp/grdblend_resampled_27780_2.nc -Vi -fg --GMT_HISTORY=false
grdsample [INFORMATION]: Processing input grid
grdsample [INFORMATION]: Round-off patrol changed geographic grid increment for latitude from 0.00416666666666666834 to 0.00416666666666666661
grdsample [INFORMATION]: Given domain implies x_inc = 0.000833333
grdsample [INFORMATION]: Given domain implies y_inc = 0.000833333
grdsample [INFORMATION]: Input  grid (0/10/40/50) n_columns = 2400 n_rows = 2400 dx = 0.00416666666667 dy = 0.00416666666667 registration = 1
grdsample [INFORMATION]: Output grid (3.96666666667/4.27083333333/44/44.25) n_columns = 366 n_rows = 301 dx = 0.000833333333315 dy = 0.000833333333333 registration = 0
grdsample [INFORMATION]: Reading grid from file ~/.gmt/server/earth/earth_relief/earth_relief_15s_p/N40E000.earth_relief_15s_p.nc
grdsample [INFORMATION]: Set boundary condition for left   edge: natural
grdsample [INFORMATION]: Set boundary condition for right  edge: natural
grdsample [INFORMATION]: Set boundary condition for top    edge: natural
grdsample [INFORMATION]: Writing grid to file /tmp/grdblend_resampled_27780_2.nc
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.000833333333315067399 to 0.000833333333333333387
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmin from 3.966666666666667 to 3.966666666666667
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmax from 4.270833333333334 to 4.270833333333334
grdblend [INFORMATION]: Blend file /tmp/grdblend_resampled_27780_2.nc in 0/10/40/50 with normal weight 1 [0-300]
grdblend [INFORMATION]: Processing input grids
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.000833333333315067399 to 0.000833333333333333387
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmin from 3.966666666666667 to 3.966666666666667
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmax from 4.270833333333334 to 4.270833333333334
grdblend [INFORMATION]: Processed row     301 of 301
grdblend [INFORMATION]: Referencing grid data to GMT_GRID memory location
grdblend [INFORMATION]: Set boundary condition for all edges: natural
grdblend [INFORMATION]: Set boundary condition for left   edge: natural
grdblend [INFORMATION]: Set boundary condition for right  edge: natural
grdblend [INFORMATION]: Set boundary condition for bottom edge: natural
grdblend [INFORMATION]: Set boundary condition for top    edge: natural
grdblend [INFORMATION]: Blended grid size of @GMTAPI@-S-O-G-G-G-Y-000002 is 366 x 301
grdblend [INFORMATION]: All nodes assigned values
grdimage [INFORMATION]: Set boundary condition for all edges: natural
grdimage [INFORMATION]: Set boundary condition for left   edge: natural
grdimage [INFORMATION]: Set boundary condition for right  edge: natural
grdimage [INFORMATION]: Set boundary condition for bottom edge: natural
grdimage [INFORMATION]: Set boundary condition for top    edge: natural
grdimage [INFORMATION]: Central meridian not given, default to 4.12
grdimage [INFORMATION]: Map scale is 1.6435 km per cm or 1:164350.
grdimage [INFORMATION]: Derive intensity grid from data grid
grdimage [INFORMATION]: Calling grdgradient with args -G@GMTAPI@-S-O-G-G-G-Y-000015 -A-45.0 -Nt1+a0 -R3.968333333333333/4.270833333333333/44/44.25 --GMT_HISTORY=false @GMTAPI@-S-I-G-G-G-N-000001
grdgradient [INFORMATION]: Processing input grid
grdgradient [INFORMATION]: Set boundary condition for right  edge: natural
grdgradient [INFORMATION]: Set boundary condition for bottom edge: natural
grdgradient [INFORMATION]: Set boundary condition for top    edge: natural
grdgradient [INFORMATION]:  Min Mean Max sigma intensities:grdgradient [INFORMATION]: 	-0.720555859751	-0.00802364892491	0.890532427816	0.1258832986
grdgradient [INFORMATION]: Referencing grid data to GMT_GRID memory location
grdgradient [INFORMATION]: Set boundary condition for right  edge: natural
grdgradient [INFORMATION]: Set boundary condition for bottom edge: natural
grdgradient [INFORMATION]: Set boundary condition for top    edge: natural
grdimage [INFORMATION]: Allocate and read data from file ~/.gmt/sessions/gmt_session.27776/=tiled_80_PO.000000
grdimage [INFORMATION]: Allocates memory and read intensity file
grdimage [ERROR]: Dimensions of intensity grid do not match that of the data grid!

The plot shows up empty, no terrain or anything:

version information:

PyGMT information:
  version: v0.1.2+40.gfc246e7a
System information:
  python: 3.7.6 | packaged by conda-forge | (default, Jun  1 2020, 18:57:50)  [GCC 7.5.0]
  executable: ~/miniconda3/envs/pygmt/bin/python
  machine: Linux-5.4.0-42-generic-x86_64-with-debian-bullseye-sid
Dependency information:
  numpy: 1.19.1
  pandas: 1.1.1
  xarray: 0.16.0
  netCDF4: 1.5.3
  packaging: 20.4
  ghostscript: 9.22
  gmt: 6.1.0
GMT library information:
  binary dir: ~/miniconda3/envs/pygmt/bin
  cores: 6
  grid layout: rows
  library path: ~/gmt-install-dir/lib/libgmt.so
  padding: 2
  plugin dir: ~/gmt-install-dir/lib/gmt/plugins
  share dir: ~/gmt-install-dir/share
  version: 6.1.1

I noticed that this is only a problem for the 3 arc second and 1 arc second earth relief grids, using 15 arc second (i.e. @earth_relief_15s) and lower resolutions don’t seem to be a problem.

pwessel · August 30, 2020, 11:57pm

Yes, the same issue was reported here: Gradient problem using earth_relief_03s
THe workaround is to give w/e/s/n in exact multiples of 3 arc seconds. The command-line script works for me with -R3:57/4:18/44:00/44:15.

weiji14 · August 31, 2020, 12:20am

Ah good, thanks Paul! I thought I had seen that problem somewhere before. @Saren, the equivalent PyGMT code would be as below (note that we’re setting the region to be in degrees, minutes, seconds instead of decimal degrees):

import pygmt

# region = [3.970, 4.270, 44.000, 44.250]
region = ["3:57", "4:18", "44:00", "44:15"]
fig = pygmt.Figure()
fig.basemap(region=region, projection="M8i", frame=True)
fig.grdimage("@earth_relief_03s", cmap="geo", shading=True)
fig.coast(rivers="a/0.5,skyblue", water="skyblue")
fig.coast(borders="1/1,black")
fig.show()

produces:

Might need to play around with the rivers/water/borders setting, doesn’t seem to show up on the map.

pwessel · August 31, 2020, 1:00am

And if you are completely on land then you can specify srtm_relief_03s instead and it will work with your original region. I am working on the underlying problem.

pwessel · August 31, 2020, 3:44am

Underlying problem has been fixed in GMT master and 6.1 branches.

Saren · August 31, 2020, 6:42am

Thanks you all for your help and fast answer!