[GIS] Merge Multi-Band Geotiff with different datatypes (with gdal)

gdalgeotiff-tiffmergemulti-band

I have three geotiffs:

  • orthomosaic.tif: a multispectral orthomosaic with three channels:

    Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
Band 4 Block=512x512 Type=Byte, ColorInterp=Undefined

  • DTM.tif (Digital Terrain Model) and DSM.tif (Digital Surface Model) of the same extent with one band each:

    Band 1 Block=5001x1 Type=UInt16, ColorInterp=Gray

I want to merge those three into one single file (to consequently retile them and work with smaller tiles having all six channels). This works fine with this command where the parameter -seperate preserves the bands and -o stands for "output":

gdal_merge.py -o merged_file.tif -separate  orthomosaic.tif  DTM.tif  DSM.tif

  • Output is: merged_file.tif

    Band 1 Block=10002x1 Type=Byte, ColorInterp=Gray
Band 2 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 3 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 4 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 5 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 6 Block=10002x1 Type=Byte, ColorInterp=Undefined

Problem is: The script converts all channels to Type=Byte, which makes the DTM/DSM channels Type=UInt16 useless for me.

Question is: Is it possible to have different data types in different bands in a single geotiff file? And how do I do this when merging several geotiffs?

Best Answer

You can try to use gdalbuildvrt (http://www.gdal.org/gdalbuildvrt.html) instead of gdal_merge.py. As the result the xml-like file describing virtual raster created. GDAL and QGIS works with this file. 

gdalbuildvrt qqv.vrt -separate qq1.tif qq3.tif

The only problem is, that only first band from all files will be added. This simply fixed by editing vrt file.

$ gdalinfo qqv.vrt 
Driver: VRT/Virtual Raster
Files: qqv.vrt
   /home/bishop/tmp/qq1.tif
   /home/bishop/tmp/qq3.tif
Size is 735, 547
Coordinate System is:
GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",
  DATUM["Not_specified_based_on_WGS_84_ellipsoid",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6030"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4030"]]
Origin = (513165.703783421427943,5016826.563432835973799)
Pixel Size = (59.969933184855272,-60.001568937759295)
Corner Coordinates:
Upper Left  (  513165.704, 5016826.563) (Invalid angle,Invalid angle)
Lower Left  (  513165.704, 4984005.705) (Invalid angle,Invalid angle)
Upper Right (  557243.605, 5016826.563) (Invalid angle,Invalid angle)
Lower Right (  557243.605, 4984005.705) (Invalid angle,Invalid angle)
Center      (  535204.654, 5000416.134) (Invalid angle,Invalid angle)
Band 1 Block=128x128 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000 
  NoData Value=nan
Band 2 Block=128x128 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000 
  NoData Value=nan
Band 3 Block=128x128 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000 
  NoData Value=nan
Band 4 Block=128x128 Type=Int16, ColorInterp=Undefined

Raster VRT format tutorial: http://www.gdal.org/gdal_vrttut.html

Related Question