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?

2 Answers 2


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, osr
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

#min & max random values of the output raster

## 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 upper left corner of the image
#in same units as spatial reference

#Cellsize in same units as spatial reference

#dst_ds.SetGeoTransform( [ xmin, cellsize, 0, ymax, 0, -cellsize ] )
#srs = osr.SpatialReference()
#dst_ds.SetProjection( srs.ExportToWkt() )

raster = numpy.random.randint(zmin,zmax, (nbands, nrows, ncols)).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
  • Thanks a lot, where can read what these things do? SetUTM (ok I know what that does) SetWellKnown GeogCS, se projection, set geo transform, etc... but looks like exactly what I need. Thanks a lot!
    – JEquihua
    Jul 11, 2012 at 16:20
  • For more info on the georeferencing parts of the code, see the Projections Tutorial - gdal.org/ogr/osr_tutorial.html
    – user2856
    Jul 11, 2012 at 23:59

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.


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

  • No, I really need to make my own image. The page is interesting, thank you, what I would really need is the missing example, how to save a 3d numpy array as a multi-band geoTIFF. But thanks a lot!
    – JEquihua
    Jul 11, 2012 at 6:53
  • Edited with more info Jul 11, 2012 at 7:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.