gdal_merge.py is the correct tool to 'stack' your input images.
Assuming that your first band has a valid color table you could use:
gdal_merge.py -separate -pct -o output_file.tif file1.tif file2.tif file3.tif
Note: The command has been reformatted with -o output_file.tif
before the list of inputs.
From the docs:
-pct:
Grab a pseudocolor table from the first input image, and use it for the output. Merging pseudocolored images this way assumes that all input files use the same color table.
I would test your output with gdalinfo -stats
to make sure it is being stacked properly.
Updated for OP
From the osgeo list, it looks like you might try a different format to check the results:
There's no way in the TIFF format to really specify the color interpretation of each band. The way GDAL builds the color interpretation when reading a TIFF file is a combination of the value of the PHOTOMETRIC and EXTRASAMPLES tag.
-Evan (the poster) knows GDAL inside and out.
i've found a way using a gdal object supported by numpy arrays..
import gdal
import numpy as np
from affine import Affine
def save_multiband (output_name, dataset, raster_data, driver,
NaN_Value, nband, arr_projection=None):
if arr_projection is None:
arr_projection = []
else:
if str(type(arr_projection)) == "<class 'osgeo.gdal.Dataset'>":
srs = arr_projection.GetProjectionRef()
elif str(type(arr_projection)) == "<type 'int'>":
srs = osr.SpatialReference()
srs.ImportFromEPSG(arr_projection)
srs = srs.ExportToPrettyWkt()
else:
srs = arr_projection
NaN_rast = NaN_Value
# Open the reference dataset
g = dataset
# Get the Geotransform vector
if raster_data is False:
raster_data = g.ReadAsArray()
rast_arr = np.copy(raster_data)
#auxiliar
band = 0
if type(raster_data) == tuple:
rast_arr = np.array(raster_data[band])
if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
geo_transform = g.GetGeoTransform()
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize # Raster ysize
srs = g.GetProjectionRef() # Projection
elif str(type(g)) == "<class 'affine.Affine'>":
geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])
rast_arr = raster_data[band,:,:]
x_size = int(rast_arr.shape[1])
y_size = int(rast_arr.shape[0])
driver = gdal.GetDriverByName(driver)
dataset_out = driver.Create(output_name, x_size, y_size, nband,
gdal.GDT_Float32)
#end auxiliar
for band in range(nband):
if type(raster_data) == tuple:
rast_arr = np.array(raster_data[band])
if str(type(g)) == "<class 'osgeo.gdal.Dataset'>":
geo_transform = g.GetGeoTransform()
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize # Raster ysize
srs = g.GetProjectionRef() # Projection
elif str(type(g)) == "<class 'affine.Affine'>":
geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])
rast_arr = raster_data[band,:,:]
x_size = int(rast_arr.shape[1])
y_size = int(rast_arr.shape[0])
#PROCESS RASTERIO NUMPY
else:
geo_transform = (g[1][2], g[1][0], g[1][1], g[1][5], g[1][3], g[1][4])
rast_arr = np.array(g[0])
x_size = int(rast_arr.shape[2])
y_size = int(rast_arr.shape[1])
rast_arr[rast_arr == NaN_rast] = np.NaN
dataset_out.SetGeoTransform(geo_transform)
dataset_out.SetProjection(srs)
dataset_out.GetRasterBand(band + 1).WriteArray(
rast_arr.astype(np.float32))
using a example data:
#assuming the data had the same extent and spatial resolution
B1 = gdal.Open(r'C\B1.tif')
B2 = gdal.Open(r'C\B2.tif')
B3 = gdal.Open(r'C\B3.tif')
B4 = gdal.Open(r'C\B4.tif')
srs = B1.GetProjectionRef()
b1, b2, b3, b4 = (B1.ReadAsArray(), B2.ReadAsArray(),
B3.ReadAsArray(), B4.ReadAsArray())
def GetTransform(Transform):
transform = Transform.GetGeoTransform()
Aff = Affine(transform[1], transform[2], transform[0], transform[4],
transform[5], transform[3])
return(Aff)
B = np.stack([b1, b2, b3, b4])
save_multiband(path+'\\SR.tif', GetTransform(B1), B, "GTiff", -9999, B.ndim + 1, srs)
Best Answer
You can use a simple python script that uses the
Band.SetDescription
method to set the band names:This functionality should probably be in
gdal_edit
, but isn't.For example: