Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif.
From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it in QGIS. If I remember correctly, when you first add a raster, it often shows up all one color if you don't have a pre-existing color map.
import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
print 'Could not open image file'
sys.exit(1)
# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize
cropData = band1.ReadAsArray(0,0,cols,rows)
listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]
# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
print 'Could not create reclass_40.tif'
sys.exit(1)
outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)
for i in range(0, rows):
for j in range(0, cols):
if cropData[i,j] in listAg:
outData[i,j] = 100
elif cropData[i,j] in listNotAg:
outData[i,j] = -100
else:
outData[i,j] = 0
# write the data
outBand.WriteArray(outData, 0, 0)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
del outData
I'm not totally certain on what you need, but I think:
gdalbuildvrt -separate stack.vrt lsat1.tif lsat2.tif ...
Should give you a gdal dataset that is 'stacked and is covers the extent of all images. If you need a tif after that, use
gdal_translate stack.vrt stack.tif
Best Answer
May be this is what you are looking for:
Using opencv library
img_p is the panchromatic image img_ms is the multispectral image This will resize the multispectral image to the exact dimensions of the panchromatic.