Exporting embellished rasters with proper coordinates

Hi all,

Condensed problem: I want to export a georeferenced geotiff raster with contours.

I have a bathymetry grid (geotiff). The coordinates are in projected units.

If I use

gmt grdimage -Abat_gmt.tif bat.tif -I+d

I get a nice color-relief map with shading. The coordinates are preserved and I can plot it in e.g. QGIS.

  • NaN values are represented with the value 99999 in my grid. How can I tell grdimage to use this as a NaN-value?

I would like to plot contours as well. That means I need to bake a GMT cake, calling e.g. grdcontour and other modules (as needed). After applying contours I want to export it, keeping the proper coordinates, like the result in the grdimage call above.

Any tips on how doing this?

I tried the psconvert approach (https://docs.generic-mapping-tools.org/latest/psconvert.html#w), adding a %%PROJ line in the .ps-file to specify the projection that should be used. But I’m, confused what to use when the grid is in already in projected units.

To test with geodetic coordinates, I’ve inverse projected the grid, and did a grdimage, and then injected:

%%PROJ utm 6.663782868 8.589722967 73.280813537 73.680863706 238337 295263 8149712 8187663 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

in the first line of test.ps and ran gmt psconvert -Tt -W+g test.ps (this creates two files; test.tif and test.tiff - should it do that?)

$ gmt psconvert -Tt -W+g test.ps -V
psconvert [INFORMATION]: Processing test.ps...
psconvert [INFORMATION]: Find HiResBoundingBox ...
psconvert [INFORMATION]: Figure dimensions: Width: 443.772 points [15.6553 cm]  Height: 607.74 points [21.4397 cm]
psconvert [INFORMATION]: [79.236 72 523.008 679.74]...
psconvert [INFORMATION]: Convert to TIF...
psconvert [INFORMATION]: width = 2533	height = 1850	X res = 22.473747	Y res = 20.525149
psconvert [INFORMATION]: Wrote world file test.tfw
psconvert [INFORMATION]: Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Input file size is 2533, 1850
0...10...20...30...40...50...60...70...80...90...100 - done.
psconvert [INFORMATION]: The gdal_translate command:
gdal_translate -mo TIFFTAG_XRESOLUTION=300 -mo TIFFTAG_YRESOLUTION=300 -a_srs '+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs' -co COMPRESS=LZW -co TILED=YES  'test.tif' 'test.tiff'

The resulting tiff has reasonable y-coordinates, but the x-coordinates are way off:

$ gmt grdinfo test.tif
test.tif: Title: Grid imported via GDAL
test.tif: Command:
test.tif: Remark:
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: -96523725866.8 x_max: -96523668940.8 x_inc: 22.4737465456 name: x n_columns: 2533
test.tif: y_min: 8149712 y_max: 8187683.52515 y_inc: 20.525148729 name: y n_rows: 1850
test.tif: v_min: 24 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0

= Other approaches

== mapnik
I read some blog posts about mapnik; seemed to be a nice solution. You create an .xml-file with all the layers you want (tif’s, shapefiles etc) and mapnik renders a final nice image with all layers combined.

I tried this, but setting up mapnik, for me, was a terrible experience. All kind of python and gdal plugins needing to be manually compiled. It made a mess of my system, with no result.

== gdal
I tried using a workflow based on gdal; 1)make color relief, 2)make hillshade, 3)combine these, 4)extract contours and 5)burn contours into the final raster.

This gives an OK result, but the contours becomes very rough and blocky. This is related to that the finest features cannot be thinner/finer than the raster resolution. One solution is to resample the grid to an artificial higher resolution, so that the contours become nicer, but this cause huge rasters.

gdal has many features that I’ve neved (had to) used, but they are quite handy to have in the toolbox. See below for the gdal workflow I used.

#make cpt
gmt makecpt -Cabyss -T-4000/-2000/50 -N | awk '{print $1, $2}' | sed 's/\// /g' > col.cpt

#resample grid if necessary

#make a color-relief dem based on the cpt created above
gdaldem color-relief -b 1 bat.tif col.cpt bat_col.tif

#make shade
gdaldem hillshade bat.tif bat_shade.tif

#gamma hillshade
gdal_calc.py -A bat_shade.tif --outfile=bat_shade_gamma.tif --calc="uint8(((A / 255.)**(1/0.5)) * 255)" --overwrite

# overlay shade and color-relieff
gdal_calc.py -A bat_shade_gamma.tif -B bat_col.tif --allBands=B \
--calc="uint8( ( \
                 2 * (A/255.)*(B/255.)*(A<128) + \
                 ( 1 - 2 * (1-(A/255.))*(1-(B/255.)) ) * (A>=128) \
               ) * 255 )" --outfile=bat_final.tif --overwrite

#extract contours
gmt grdcontour bat.tif -D%c.gmt -C5 -Q20

#add OGR-header. gdal requires this. @Je is the epsd code of the projection
sed -i '1i # @VGMT1.0 @GLINESTRING # @Je25833' C.gmt

#burn in the contours to the raster
gdal_rasterize -b 1 -b 2 -b 3 -burn 0 -burn 0 -burn 0 C.gmt bat_final.tif

#export to tiles, if wanted
#gdal2tiles.py bat_final.tif

Andreas,

You should not have to mess around with the %%PROJ comment line in the postscript. That line is generated by GMT. You should only have to use psconvert -W+g ... and get the geotiff. Doesn’t it work?

Geotiff files have a nodata variable whose content is detected by GMT and used instead of the default NaN. Also here, didn’t it work?

Joaquim,

Thank you for your interest.

I’ll go through the steps.

My starting point is my well behaved tif-file. With valid coordinates;

$ gmt grdinfo BathymetryExport-dtm.tif 
BathymetryExport-dtm.tif: Title: Grid imported via GDAL
BathymetryExport-dtm.tif: Command: 
BathymetryExport-dtm.tif: Remark: 
BathymetryExport-dtm.tif: Pixel node registration used [Cartesian grid]
BathymetryExport-dtm.tif: Grid file format: gd = Import/export through GDAL
BathymetryExport-dtm.tif: x_min: 238337.5 x_max: 295262.5 x_inc: 25 name: x n_columns: 2277
BathymetryExport-dtm.tif: y_min: 8149712.5 y_max: 8187662.5 y_inc: 25 name: y n_rows: 1518
BathymetryExport-dtm.tif: v_min: -3303.82006836 v_max: -1646.84997559 name: z
BathymetryExport-dtm.tif: scale_factor: 1 add_offset: 0
+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Doing a gmt grdimage BathymetryExport-dtm.tif -Agrdimage.tif

Gives me a color-relief and correct georeferencing:

$ gmt grdinfo grdimage.tif 
grdimage.tif: Title: Grid imported via GDAL
grdimage.tif: Command: 
grdimage.tif: Remark: 
grdimage.tif: Pixel node registration used [Cartesian grid]
grdimage.tif: Grid file format: gd = Import/export through GDAL
grdimage.tif: x_min: 238337.5 x_max: 295262.5 x_inc: 25 name: x n_columns: 2277
grdimage.tif: y_min: 8149712.5 y_max: 8187662.5 y_inc: 25 name: y n_rows: 1518
grdimage.tif: v_min: 24 v_max: 254 name: z
grdimage.tif: scale_factor: 1 add_offset: 0
+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

NaN’s in the tif turns into RGB 128/128/128 (gray) in the tif.

If I run grdimage with -Q (“Make grid nodes with z = NaN transparent”), I still get 128/128/128 in the area with NaN’s, but the tif gets a fourth band;

$ gdalinfo BathymetryExport-dtm.tif  | grep NoData
  NoData Value=99999

gdalinfo grdimage.tif | grep NoData
  NoData Value=0
  NoData Value=0
  NoData Value=0

$ gdalinfo grdimage-with-Q.tif | grep NoData
  NoData Value=0
  NoData Value=0
  NoData Value=0
  NoData Value=0

So the original tif has one NoData value (applying for all bands?) while the grdimage generated ones have one NoData value for each band (including the fourth band with -Q) (but it’s 0 for all of them)

Trying the simple script

gmt grdimage -JX10c BathymetryExport-dtm.tif > test.ps
gmt psconvert -Tt+w test.ps

Creates this tiff:

 gmt grdinfo test.tif 
test.tif: Title: Grid imported via GDAL
test.tif: Command: 
test.tif: Remark: 
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: 0 x_max: 2480 x_inc: 1 name: x n_columns: 2480
test.tif: y_min: 0 y_max: 3509 y_inc: 1 name: y n_rows: 3509
test.tif: v_min: 24 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0

Note coordinates - not in the original system anymore. Should this be working?

Which means it working fine. You got the default NaN color of the CPTs. For other color you have to create your own cpt with a different NaN color.

Hmm, what it -Tt+w? Isn’t -Tt -W+g what you want?

I’m sorry. Yes, of course it is.

gmt grdimage -JX10c BathymetryExport-dtm.tif > test.ps
gmt psconvert -Tt -W+g -V test.ps
psconvert [INFORMATION]: Processing test.ps...
psconvert [INFORMATION]: Find HiResBoundingBox ...
psconvert [INFORMATION]: Figure dimensions: Width: 283.448 points [9.99942 cm]  Height: 283.446 points [9.99934 cm]
psconvert [INFORMATION]: [239.56 72 523.008 355.446]...
psconvert [ERROR]: Requested an automatic geotiff generation, but no recognized psconvert option was used for the PS creation.
psconvert [INFORMATION]: Convert to TIF...
psconvert [INFORMATION]: width = 1182	height = 1182	X res = 48.159898	Y res = 32.106599
psconvert [INFORMATION]: Wrote world file test.tfw
psconvert [ERROR]: Could not find the Proj4 command in the PS file. No conversion performed.

That’s why I injected proj-stuff in my .ps. Should it work without having to do so?

There is a file created, test.tif:

gmt grdinfo test.tif
test.tif: Title: Grid imported via GDAL
test.tif: Command: 
test.tif: Remark: 
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: -206844860108 x_max: -206844803183 x_inc: 48.1598984772 name: x n_columns: 1182
test.tif: y_min: 8149712.5 y_max: 8187662.5 y_inc: 32.1065989848 name: y n_rows: 1182
test.tif: v_min: 24 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0

The y-coordines looks OK, but the x-coordinates are not.

OK, I can also confirm the

psconvert [ERROR]: Requested an automatic geotiff generation, but no recognized psconvert option was used for the PS creation.
psconvert [ERROR]: Could not find the Proj4 command in the PS file. No conversion performed.

it must be some psconvert regression, but despite that it still worked in my example. What I did was to extract a bit of the central Atlantic, convert it to UTM and save in geotif.

gmt grdimage plateau.tiff -JX15 -P -K > p.ps
gmt grdcontour plateau.tiff -C1000 -J -O >> p.ps
psconvert p.ps -Tt -W+g

Next I assigned it the proj info +proj=utm +zone=26 +datum=WGS84 +units=m +no_defs in Mirone and plotted the coastlines. The match was as good as one can visually tell. So, my conclusion is that it still works.

Joaquim,

I tried your commands with my tif, and it works!

It seems that the choice of -J is important.

gmt grdimage BathymetryExport-dtm.tif -Jx1:10000 -P --PS_MEDIA=A0 -I+d > p.ps
gmt psconvert p.ps -Tt -W+g -V
gmt grdinfo p.tif
p.tif: Title: Grid imported via GDAL
p.tif: Command: 
p.tif: Remark: 
p.tif: Pixel node registration used [Cartesian grid]
p.tif: Grid file format: gd = Import/export through GDAL
p.tif: x_min: 238336.653343 x_max: 246478.956087 x_inc: 0.84665724697 name: x n_columns: 9617
p.tif: y_min: 3644541941.56 y_max: 3644553569.64 y_inc: 0.84666354326 name: y n_rows: 13734
p.tif: v_min: 1 v_max: 255 name: z
p.tif: scale_factor: 1 add_offset: 0

The y-coordinates here are wrong.

Change the scale (add another 0), and it works:

gmt grdimage BathymetryExport-dtm.tif -Jx1:100000 -P --PS_MEDIA=A0 -I+d > p.ps

gmt psconvert p.ps -Tt -W+g -V

gmt grdinfo p.tif
p.tif: Title: Grid imported via GDAL
p.tif: Command: 
p.tif: Remark: 
p.tif: Pixel node registration used [Cartesian grid]
p.tif: Grid file format: gd = Import/export through GDAL
p.tif: x_min: 238329.034057 x_max: 295254.034057 x_inc: 8.46594289114 name: x n_columns: 6724
p.tif: y_min: 8149712.5 y_max: 8187662.5 y_inc: 8.4653134062 name: y n_rows: 4483
p.tif: v_min: 1 v_max: 255 name: z
p.tif: scale_factor: 1 add_offset: 0

Why would this make a difference?

Don’t know. Could it be that the larger scale map overflows the paper size and gets clipped and thus completely fooling the algorithm that converts points to meters?

I think you are correct. If the image is not able to fit the (chosen) PS_MEDIA, the coordinates (y-coordinates in my case) blows up.

@Andreas, any more info where that gdal_calc equation is coming from? The condition on < > 128 makes me think on a illuminating/shading condition but the rest is a bit fuzzy for me.

@Joaquim, I found it while googling how to merge color-relief and hillshade. See https://gis.stackexchange.com/questions/255537/merging-hillshade-dem-data-into-color-relief-single-geotiff-with-qgis-and-gdal/255574#255574.

Thanks. Will try when time permits. Should be a one-liner in Julia.

Another one on this.

-P seems to be important:

This grid is correctly placed (with -P):

gmt grdimage BathymetryExport-dtm.tif -JX10c -P --PS_MEDIA=A0 > test.ps
gmt psconvert -Tt -W+g test.ps

gmt grdinfo test.tif
test.tif: Title: Grid imported via GDAL
test.tif: Command: 
test.tif: Remark: 
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: 238289.340102 x_max: 295214.340102 x_inc: 48.1598984772 name: x n_columns: 1182
test.tif: y_min: 8149712.5 y_max: 8187662.5 y_inc: 32.1065989848 name: y n_rows: 1182
test.tif: v_min: 255 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0

Omit the -P and the x-coordinates are destroyed:

gmt grdimage BathymetryExport-dtm.tif -JX10c --PS_MEDIA=A0 > test.ps
gmt psconvert -Tt -W+g test.ps

gmt grdinfo test.tif
test.tif: Title: Grid imported via GDAL
test.tif: Command: 
test.tif: Remark: 
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: -206844501943 x_max: -206844445018 x_inc: 48.1598984772 name: x n_columns: 1182
test.tif: y_min: 8149744.6066 y_max: 8187662.5 y_inc: 32.1065989848 name: y n_rows: 1181
test.tif: v_min: 255 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0