How to plot DEM grid using -JG (perspective) projection?

Q&A
1

Hi!

I am having some problems trying to use grdimage (or plotting a DEM) using -JG.

In past versions of GMT I have used the “-p” option in grdimage to make perspective plots with an overlain grid on my basemap.

However, this time, in GMT 6.5, I am making a plot that is similar-ish to
this example (bottom one); so using the perspective projection “-JG” for a large slice of the planet.

Now, my script runs fine when I am just using ‘coast’ like the example shows

PROJ=-JG${longitude}/${latitude}/12c+z${altitude}+a${azimuth}+t${tilt}+w${twist}+v${Width}/${Height}
gmt basemap $PROJ -Rg -BESnw+gslategray -By10g5 -Ba20xg10
gmt coast $PROJ -By10g5 -Ba20xg10 -Glightbrown -Slightblue -W -Ia/blue -Di -Na -Rg

But when I try to switch ‘coast’ to a grid being displayed (specifically the DEM of Mars which is what need to be plotting):

gmt grdimage ./grids/mars_dem.tif -C’matlab/copper’ -t35 -I+a-45+nt0.5+m0

I get the following error:

grdimage [ERROR]: gmt_gdalread: Computed -srcwin falls outside raster size of 720x360 in the x-direction.
grdimage [ERROR]: ERROR reading file with gdalread.
grdimage (gmtapi_import_grid): Could not open file [./grids/mars_dem.tif]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

By the way, when I add some points using gmt plot, it also works perfectly fine, it’s just the grdimage that gmt doesn’t like in this perspective mode.

The DEM/grid file is fine by the way, I plot it in the same subfigure just before for a global map using “-JN”, so the problem is the “-JG” projection or something else.

Happy to provide more detail/code/variables/images/data if and when needed.

Many thanks for any help in advance!

2

We need a reproducible command (there is no -p in the command that you say it fails). And no $$${…} please. For example, I use Windows.

3

Sure! Here is a reduced form of the script:

#!/bin/bash

# Set parameters for the most beautiful plot
gmt set PS_LINE_CAP = round
gmt set PS_LINE_JOIN = miter
gmt set FONT_LABEL 16p,Helvetica,Black
gmt set FONT_ANNOT_PRIMARY 12p,Helvetica,Black
gmt set FORMAT_GEO_OUT = D
gmt set FORMAT_GEO_MAP = D
gmt set PS_MEDIA = a2
gmt set FONT_TAG = auto,Helvetica-Bold
gmt set MAP_FRAME_TYPE plain
gmt set PS_PAGE_ORIENTATION = landscape
gmt set PROJ_ELLIPSOID = Mars

#plotting options:
latitude=-80; longitude=0; altitude=1000; tilt=45; azimuth=210; twist=60; Width=20; Height=20
mapheight=3i
gmt begin mars_test
    gmt subplot begin 1x2 -Fs5i,5i/5i -M0.2i -Ad+gwhite+v+jTL -Ba0
        gmt subplot set 
        PROJ1="-JN0/${mapheight}+dh"
        gmt basemap $PROJ1 -Bnswe+gslategray -Rg
        gmt grdimage $PROJ1 ./mars_dem.tif -C'matlab/copper' -t35 -I+a-45+nt0.5+m0

        gmt subplot set
        PROJ2=-JG${longitude}/${latitude}/12c+z${altitude}+a${azimuth}+t${tilt}+w${twist}+v${Width}/${Height}
        gmt basemap $PROJ2 -Rg -BESnw+gslategray -By10g5 -Ba20xg10
        #gmt grdimage $PROJ2 ./mars_dem.tif -C'matlab/copper' -t35 -I+a-45+nt0.5+m0 -BESnw #doesnt work

    gmt subplot end
gmt end show

I’ve uploaded the DEM as well!
mars_dem.tif (507.2 KB)

I can plot the first subfigure with JN, but the second one doesnt work with JG… I’ve tried grdview, but if I have to supply a -p I am not sure how the variables would be set to match?

4

No variables (${...}) and avoid bash scripts remember? One-liners are enough.

gmt grdimage  mars_dem.tif -t35 -I+a-45+nt0.5+m0 -BESnw -JG0/-80/12c+z1000+a210+t45+w60+v20/20 -png lixo

GDAL is not able to recognize that projection. Convert the tif file to netcdf with grdconvert and it will work.

5

Adding:

gmt grdconvert ./grids/mars_dem.tif -G./grids/mars_dem.nc

now makes it work! Muito obrigada!