Creating a multispectral image from scratch

Question

I want to make a multispectral image from cero to do some tests on it. Something really simple like 5 completely uniform bands with salt and pepper noise on them or a square of different values at the center. Clearly this would just be a stack of matrices, a multidimensional array, which is pretty straight forward to generate. I want to achieve this using python and gdal but gdal is being pretty hermetic, I don't get the hang of it at all. I ideally would want to create a geotiff file. Could anyone help me on this? some pointers or some gdal tutorial which is very gentle? Thank you all.

Answer 1

You want the gdal.band.WriteArray method. There's an example in the GDAL API tutorial (reproduced below):

format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )
dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )

srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
srs.SetWellKnownGeogCS( 'NAD27' )
dst_ds.SetProjection( srs.ExportToWkt() )

raster = numpy.zeros( (512, 512), dtype=numpy.uint8 )    
dst_ds.GetRasterBand(1).WriteArray( raster )

# Once we're done, close properly the dataset
dst_ds = None

For generating the random data,look at the numpy.random module.

Here's a more complete working example:

from osgeo import gdal
import numpy

dst_filename = '/tmp/test.tif'
#output to special GDAL "in memory" (/vsimem) path just for testing
#dst_filename = '/vsimem/test.tif'

#Raster size
nrows=512
ncols=512
nbands=7

#min & max random values of the output raster
zmin=0
zmax=12345

## See http://gdal.org/python/osgeo.gdal_array-module.html#codes
## for mapping between gdal and numpy data types
gdal_datatype = gdal.GDT_UInt16
np_datatype = numpy.uint16

driver = gdal.GetDriverByName( "GTiff" )
dst_ds = driver.Create( dst_filename, ncols, nrows, nbands, gdal_datatype )

## These are only required if you wish to georeference (http://en.wikipedia.org/wiki/Georeference)
## your output geotiff, you need to know what values to input, don't just use the ones below
#Coordinates of the lower left corner of the image
#in same units as spatial reference
#xllcorner=147.2  
#yllcorner=-34.54

#Cellsize in same units as spatial reference
#cellsize=0.01

#dst_ds.SetGeoTransform( [ xllcorner, cellsize, 0, yllcorner, 0, -cellsize ] )
#srs = osr.SpatialReference()
#srs.SetWellKnownGeogCS("WGS84")
#dst_ds.SetProjection( srs.ExportToWkt() )

raster = numpy.random.randint(zmin,zmax, (ncols, nrows,nbands)).astype(np_datatype )  
for band in range(nbands):
    dst_ds.GetRasterBand(band+1).WriteArray( raster[:,:,band] )

# Once we're done, close properly the dataset
dst_ds = None

Answer 2

I know it's not what you asked for, but if all you want is multispectral or hyperspectral sample data - this test data for the Opticks project might work. Alternately, you can get LANDSAT data directly from Earth Explorer.

This site has example code to convert a 2D numpy array to a single-band geoTIFF, and a multi-band geoTIFF to a 3D numpy array.

EDIT:

Further research finds a page of example code with the 'missing example', 3D numpy array -> multi-band geoTIFF.