How to plot satellite optical image of Geotif format and show the true color?

Hello. I have a question on the Geotif image plotting.

I hope to plot a figure using the satellite remote sensing image (Geotif format) as a background. Specially,I would like to show the real color of the optical image. I know the grdimage can read Geotif directly but do not know how to set the out image with true or false color (RGB merged or one band only). The psimage can read raster images and plot color figures, but seems can not cooperate with the map projection well (like the figures in earlier post).

So could someone give me a suggestion to do this job?

Searching the older GMT forums, I found a similar post ( But seems the grdimage has been improved now and the answer is not suitable for present.

Here is the code:

#!/usr/bin/env bash

data2=../output.png # output from gdal_translate

gmt pscoast -R$data -J$J -Df -W -K -Ba+t  > $ps
grdimage  $data -R -J -K -O >> $ps

gmt psimage $data2 -Dx0/0+w10c -J -O -X12c  -Ba2 -R >> $ps
gmt psconvert $ps -A -P -Tg

Here is the figure

I just found the method to plot my figure.

grdimage ../red.tif ../green.tif ../blue.tif -R -J$J -O >> $ps

But firstly I need to split the image into three individual band files. Maybe there should have a easy way just using one single multi-band file?

Have you looked at the CookBook’s 3.26.1. Reading multi-band images?

Thanks @Joaquim.
The cookbook is so huge. And I missed that. (The grdimage man page dose not mention the related information? ).

I tried to use the +b option.

# data has four bands
gmt grdimage $data=gd+b1 -R$data -J$J -O >> $ps

But grdimage gives the message:

grdimage [ERROR]: Using data type other than byte (unsigned char) is not implemented
grdimage (api_import_image): Wrong flag passed to gmt_dist_array [../../../area_cal_low.tif]
[Session gmt.exe (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)
[Session gmt.exe (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

However, the +b in grdmath , grdinfo works well.

I use the GMT Version 6.0.0 [64-bit] under MINGW64,Win10.

I did not find where to upload a file in this website. So data was uploaded to other site. Get it Here

Sorry, can’t do anything with that link.
Regarding the message. So it probably means that your layers are unit16 (2 bytes). If that’s the case then they first need to be converted to 1 byte images and those can finally be stacked in a RGB order and produce an image. GMT doesn’t have a tool to do this directly. You will need to use grdmath. Using the Julia wrapper and some scripting would probably easier.
Even easier would be to use Mirone for that task (but needs update after installing).

Sorry, Please try this link for data and this link for the script. The tif file is re-sampled to a low resolution, so the size is about 1M. I hope they are usable.

You are right. I translate the Geotif to Byte type using GDAL (simple command) . Then the grdimage message turns to:

./ line 18:  2028 Segmentation fault  gmt grdimage $data=gd+b1 -R$data -J$J -O >> $ps

Still some problem. But the earlier error related the data type has been resolved. Is the victory closer?

Victory is always guarantied here.

I don’t get a crash but this message instead and that I’ll have to investigate. And note that I’m using GMTdev instead of 6.0

$ gmt grdimage area_cal_low_byte.tif=gd+b1 -Jm120:30/36/1:200000 >
grdimage [WARNING]: The image memory layout (BRBa) is of a wrong type. It should be BRPa.

but this works fine (except that I don’t know the order of the bands to produce a true RGB)

gmt grdimage area_cal_low_byte.tif+b2,1,0 -JM14 -Ba -png map

Really appreciate!
The band order is right and you give the true RGB figure. I just guess that the crash I meet is related to the GDAL lib or programs included in GMT*.exe Installer under Windows operation system? I will also try to run it under Linux when back to office.

I remember that I did fix something in the +b option sometime ago but it was something related to usage of bands > 9.

I use Windows too. Going to unix wont make this work for you. My working command should work for you as well and you can run it from the Win cmd shell.

Hi, GMT people

I find that the new version of GMT6.1.1 dose not work well to plot RGB of Geotiff image. So there maybe bugs about this function in GMT 6.1.1.

The code is:

#!/usr/bin/env bash
# gmt grdinfo $data+b0
gmt grdimage $data+b2,1,0 -R$data -J$J -K >$ps
gmt pscoast -R$data -J$J -Df -W  -Ba  -O >> $ps
gmt psconvert $ps -A -P -Tf

Using GMT 6.1.1 I got the wrong figure:

IF use the GMT 6.1.0, I got the right figure.

Hi, sorry for delay but we also need the data to reproduce the issue.

Hi, Joaquim
The data link is

area_cal_low.tif is a grid, not an image so you can’t select bands like that. With GMTdev this is what I get.

grdimage area_cal_low.tif+b2,1,0 -Jm120:30/36/1:200000 -png map
grdimage [ERROR]: Using this data type (Float32) is not implemented
grdimage (gmtapi_import_image): Bad measurement unit.  Choose among c|i|p [area_cal_low.tif]
Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

This is caused by the byte setting. It can be corrected by changing Float32 to Byte. Then you can use GMT properly to plot geoFIFF without such errors.
However, my finding is that GMT6.1.1 and GMT6.1.0 got different plots. And GMT 6.1.0 is correct.

#!/usr/bin/env bash

# change to Byte from Float32
gdal_translate -ot Byte -q area_cal_low.tif area_cal_low_byte.tif

# grdinfo works well on -b
gmt grdinfo $data+b0

# output png
gmt grdimage area_cal_low_byte.tif+b2,1,0 -JM14 -Ba -png map

Yes, some regression. Even if I specify +b0 it still reads 4 bands and looks like it’s using the 4rth band as a transparency.

Fixed in master. Thanks for pointing it out.

I will try later.