"ERROR - Unable to create job file" when using grdview

I’m trying to drape an georeferenced image on a map with grdview:

gmt grdview @earth_relief_01s \
  -JL-24.9/16.55/16.3/16.7/20c -R-25.14/16.75/-24.8/16.95r \
  -Ba10mg10m -I+ -Gcaboverde.tif -pdf

This fails and I only get:

gmt [ERROR]: ERROR - Unable to create job file 
/Users/kristof/.gmt/sessions/gmt6.10220/=srtm1OP.000000

You can get caboverde.tif at my server (~12 MB) for your experimentation.

What am I doing wrong here?

Thank you and all the best,
Kristof

Bug maybe. Don’t remember what docs say but you need to provide a file name after the -pdf

Hi @Joaquim this solves my initial problem – my mistake obviously. Thank you for pointing it out. But now I run into another problem. The execution just stops without returning to a prompt. Maybe another PEBKAC by me?

Full script:

#!/usr/bin/env bash

function gmt() { /usr/local/opt/gmt_dev/bin/gmt "$@"; }
export -f gmt

echo "Using GMT version: "
gmt --version

gmt grdview @earth_relief_01s \
  -JL-24.9/16.55/16.3/16.7/20c -R-25.14/16.75/-24.8/16.95r \
  -Ba10mg10m -I+ -Gcaboverde.tif -pdf draped_img

And the result:

mbp15-kristof:drape_img kristof$ ./drape_img.sh
Using GMT version:
6.1.0_9c88924_2020.05.06
grdview [NOTICE]: Earth Relief at 1x1 arc seconds tiles provided by SRTMGL1 (land only) [NASA/USGS].

grdblend [WARNING]: File @earth_relief_15s has different increments (0.00416666666667/0.00416666666667) than the output grid (0.000277777777778/0.000277777777778) - must resample
grdblend [WARNING]: File @earth_relief_15s coordinates are phase-shifted w.r.t. the output grid - must resample
grdsample [INFORMATION]: Processing input grid
grdsample [INFORMATION]: Given domain implies x_inc = 0.000277778
grdsample [INFORMATION]: Given domain implies y_inc = 0.000277778
grdsample [INFORMATION]: Input  grid (-180/180/-90/90) n_columns = 86400 n_rows = 43200 dx = 0.00416666666667 dy = 0.00416666666667 registration = 1
grdsample [INFORMATION]: Output grid (-25.1416666667/-24.8/16.75/16.95) n_columns = 1231 n_rows = 721 dx = 0.000277777777805 dy = 0.000277777777778 registration = 0
grdsample [INFORMATION]: Reading grid from file earth_relief_15s.grd
grdsample [INFORMATION]: All boundaries set via extended data.
grdsample [INFORMATION]: Output grid extrema [-1659.66/639.847] exceeds extrema of input grid [-1621/638]; to clip output use -n...+c
grdsample [INFORMATION]: Writing grid to file /var/folders/tb/cbcbcd0n6ts2bmx5vz56w_3h0000gn/T/grdblend_resampled_12098_2.nc
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.00027777777780487817 to 0.000277777777777777778
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for latitude from 0.000277777777777776802 to 0.000277777777777777778
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmin from -25.14166666666667 to -25.14166666666667
grdblend [INFORMATION]: Blend file /Users/kristof/.gmt/server/earth_relief_15s.grd in -180/180/-90/90 with normal weight 1 [0-720]
grdblend [INFORMATION]: Processing input grids
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.00027777777780487817 to 0.000277777777777777778
grdblend [INFORMATION]: Round-off patrol changed geographic grid increment for latitude from 0.000277777777777776802 to 0.000277777777777777778
grdblend [INFORMATION]: Round-off patrol changed grid limit for xmin from -25.14166666666667 to -25.14166666666667
grdblend [INFORMATION]: Processed row     721 of 721
grdblend [INFORMATION]: Referencing grid data to GMT_GRID memory location
grdblend [INFORMATION]: Blended grid size of @GMTAPI@-S-O-G-G-G-Y-000002 is 1231 x 721
grdblend [INFORMATION]: All nodes assigned values
grdview [INFORMATION]: Set boundary condition for all edges: natural
grdview [INFORMATION]: Set boundary condition for left   edge: natural
grdview [INFORMATION]: Set boundary condition for right  edge: natural
grdview [INFORMATION]: Set boundary condition for bottom edge: natural
grdview [INFORMATION]: Set boundary condition for top    edge: natural
grdview [INFORMATION]: Map scale is 1.8121 km per cm or 1:181210.
grdview [INFORMATION]: Derive intensity grid from data grid
grdview [INFORMATION]: Calling grdgradient with args /Users/kristof/.gmt/sessions/gmt6.12096/=srtm1OP.000000 -G@GMTAPI@-S-O-G-G-G-Y-000015 -A-45.0 -Nt1 -R-25.14166666666667/-24.8/16.75/16.95 --GMT_HISTORY=false
grdgradient [INFORMATION]: Processing input grid
grdgradient [NOTICE]: Earth Relief at 1x1 arc seconds tiles provided by SRTMGL1 (land only) [NASA/USGS].

grdgradient [INFORMATION]: Assembling SRTM grid from 1x1 degree tiles given by listfile /Users/kristof/.gmt/sessions/gmt6.12096/=srtm1OP.000000
grdblend [INFORMATION]: Reading Data Table from Standard Input stream
grdblend [INFORMATION]: Given domain implies x_inc = 0.000277778
grdblend [INFORMATION]: Given domain implies y_inc = 0.000277778

Yep, mine is topping at same place. It seems to be waiting for some input.

I think I provided all options necessary. What kind of input am I missing in your opinion, @Joaquim ?

It stalls at grdblend so the problem is probably there. Can you try to split the command in two? One that generates the grid and next use drape the image on it?

Not quite sure what you want me to do – grdblend is not called by me.

I think it is due to the fact the SRTM1 data is only available on land and therefore it is automatically merged with the 15s dataset to have some ocean floor. I don’t know how to split this?

SUggest you specify -Qi with some dpi. You are getting the default -Qs or -Qm by default and there seems to be a bug in there (it is using the larger dimension of the image nx, ny with the smaller dimensions of the Topo->x and Topy->y. Not clear yet. Also the -I+d leads to an issue in reading the srtm list of tiles - it hangs in gets and I think I know why. So no -I for now and maybe -Qi300.

@Joaquim could you see if there is some assumption in grdview that the image must match the topo resolution? We are definitively using Topo->x array (smaller) but accessing it in a loop over the nx from the image (larger)? Surely a bug.

Too many things at same time in my head. Maybe open an issue. Meanwhile this worked relatively fast

grdcut @earth_relief_01s -R-25.14/16.75/-24.8/16.95r -GCV.grd -Vl
grdview CV.grd -Ba10mg10m -I+ -Gcaboverde.tif -Vl -pdf lix0 -Qi100

@Joaquim your example is quite fast but doesn’t produce the desired result. If you look closely the shaded terrain is at another scale than the image to be draped.

According to the docs about grdview -G

The drapegrid may be of a different resolution than the reliefgrid.

Which leads me to believe that the grid for the shading can be of different resolution then the georeferenced image I’m trying to drape over it.

Well, it tries but …

grdview [INFORMATION]: Get and store projected vertices
grdview [INFORMATION]: Resampling relief grid to drape grid resolution
grdview [INFORMATION]: Resampling illumination grid to drape grid resolution
grdview [INFORMATION]: Start rasterization

Do you want me to make it an issue on GitHub?

Yes please, and if you can make the test case as small as possible (i.e., small grid, area)

Done. Find the issue at https://github.com/GenericMappingTools/gmt/issues/3287. I was able to reduce the uncompressed data needed to 930 KB including the relief grid. I hope this is small enough for your liking.