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?
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)
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:
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 > lixo.ps
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)
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.
Thanks.
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.
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)
Hi,
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
