Is there a mechanism to import existing georefeced image such as landsat or similar to GMT. I
Yes, see the CookBook on how to import GDAL subdatasets
Dear Joaquim,
Thank you for your reply. I did download Landsat GeoTIFF image with several bands (band4.tif). Running with grdimage spits out error message below:
But before that gdal info shows that the image comes with UTM coordinates (integer btye 16)
-------------------------------------------------------------------------------------------------
slight_smile:Band 1 Block=8231x1 Type=Int16, ColorInterp=Gray
Description = band 4 surface reflectance
---------------------------------------------------------------------------------------------------
*Do you know how to convert int16 to another format GMT can interpret?*
**gmt grdimage band4.tif=gd+b1 -Rband4.tif -JM6i -O > tt2.ps**
ERROR 5: band4.tif: GDALDataset::GetRasterBand(2) - Illegal band #
ERROR 10: Pointer 'hBand' is NULL in 'GDALGetRasterColorInterpretation'.
ERROR 10: Pointer 'hBand' is NULL in 'GDALGetRasterNoDataValue'.
ERROR 10: Pointer 'hBand' is NULL in 'GDALGetRasterColorInterpretation'.
ERROR 10: Pointer 'hObject' is NULL in 'GDALGetMetadataItem'.
ERROR 10: Pointer 'hBand' is NULL in 'GDALComputeRasterMinMax'.
***grdimage [ERROR]: Using this data type (Int16) is not implemented***
*grdimage (gmtapi_import_image): Wrong flag passed to gmt_dist_array [band4.tif]*
*[Session gmt (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)*
*[Session gmt (0)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)*
/Yacob
It’s not just converting. You must scale it down to the [0 255] interval. Check in GDAL to see if they have some option to do that. Otherwise, use gmtconvert to write a netcdf grid and next play with makecpt (you’ll need to) in order to select the data range that has … data.
The most correct way is to stretch the histogram but that requires other tools. Mirone can do it and so should the Julia wrapper do but I’m encountering some (new?) shits with it.
OK, a GMT-Julia recipe
using GMT
fname = "C:\\SIG_AnaliseDadosSatelite\\madeira\\LC08_L1TP_208037_20140110_20170426_01_T1_2014-01-10\\LC08_L1TP_208037_20140110_20170426_01_T1_2014-01-10_B4.TIF";
I = gmtread(fname);
# Stretch the histogram. To know the numbers I created an histogram of I (in Mirone)
# and came out with these numbers
I2 = mat2img(I, histo_bounds=[5820 9650]);
imshow(I2, figsize=8, savefig="LS_B4.png")
