Reading UTM geotiff DEM

Q&A
1

Hi there, sorry if this ends up being super basic. I’m trying to display this DEM (downsampled from original for upload) in order to do some swath profiling of fault scarps, in a similar way to example 33.

I’m currently unable to load my DEM, and get this output. The geographic ranges look good
grdimage [INFORMATION]: File dem-DEM_aligned_30m.tif reads with GDAL driver GTiff
grdimage [INFORMATION]: Read header from file dem-DEM_aligned_30m.tif
grdimage [INFORMATION]: File dem-DEM_aligned_30m.tif reads with GDAL driver GTiff
grdimage [INFORMATION]: gmt_map_setup perimeter search region: 70.65900000010298/71.03100000019242/39.05132776178719/39.36768301479938.
grdimage [INFORMATION]: Map scale is 0.730022 km per cm or 1:73002.2.
grdimage [INFORMATION]: Allocate and read data from file dem-DEM_aligned_30m.tif
grdimage (gmtapi_import_grid): Subset x range <= 0.0
[Session gmt (0)]: Error returned from GMT API: GMT_RUNTIME_ERROR (79)
[Session gmt (0)]: Error returned from GMT API: GMT_RUNTIME_ERROR (79)
[Session gmt (0)]: Error returned from GMT API: GMT_RUNTIME_ERROR (79)
psconvert [ERROR]: Unable to decode BoundingBox file /home/bwwj/.gmt/sessions/gmt_session.7279/psconvert_7344c.bb (maybe no non-white features were plotted?)

Here’s the file: dem-DEM_aligned_30m.tif (434.7 KB)

Here’s the code:
#!/bin/bash
gmt gmtset FORMAT_GEO_MAP D:mm

echo “testing loading tif in gmt”

im=“dem-DEM_aligned_30m”
proj=“-JU42/20c”
bounds=“-R70.759/39.154/70.931/39.265+r”

# downsample
# gdalwarp -r average -tr 30 30 -multi -wo NUM_THREADS=ALL_CPUS -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co TILED=YES -co NUM_THREADS=ALL_CPUS -co COMPRESS=LZW dem-DEM_aligned.tif dem-DEM_aligned_30m.tif

gmt begin test

  • gmt makecpt -Cdem1 -T1000/3000/50 -Z*
  • gmt grdimage ${im}.tif $bounds $proj -V #-I${im}.int -Bxa0.05f0.01 -Bya0.05f0.01*
    gmt end show
2

I’m getting this info from your grid file:

$ gmt grdinfo dem-DEM_aligned_30m.tif 
dem-DEM_aligned_30m.tif: Title: Grid imported via GDAL
dem-DEM_aligned_30m.tif: Command: 
dem-DEM_aligned_30m.tif: Remark: 
dem-DEM_aligned_30m.tif: Pixel node registration used [Cartesian grid]
dem-DEM_aligned_30m.tif: Grid file format: gd = Import/export through GDAL
dem-DEM_aligned_30m.tif: x_min: 652200.75 x_max: 665580.75 x_inc: 30 name: x n_columns: 446
dem-DEM_aligned_30m.tif: y_min: 4336968.75 y_max: 4347978.75 y_inc: 30 name: y n_rows: 367
dem-DEM_aligned_30m.tif: v_min: 1539.35083008 v_max: 2945.37670898 name: z
dem-DEM_aligned_30m.tif: scale_factor: 1 add_offset: 0
dem-DEM_aligned_30m.tif: Default CPT: 
+proj=utm +zone=42 +datum=WGS84 +units=m +no_defs

This is already projected (UTM42), so plot it with -Jx or -JX.
Not sure what the problem is. I’m getting this image from this command:

$ gmt grdimage dem-DEM_aligned_30m.tif -Jx1:100000 -png test

test
test1581×1443 238 KB

1 Like
3

Thanks for your reply, Andreas. The -Jx option worked!

1 Like
4

You’re welcome!