Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question

2 Answers 2

up vote 10 down vote accepted

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
share|improve this answer
    
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 '12 at 16:20
    
For more info on the georeferencing parts of the code, see the Projections Tutorial - gdal.org/ogr/osr_tutorial.html –  Luke Jul 11 '12 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.

EDIT:

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

share|improve this answer
    
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 '12 at 6:53
    
Edited with more info –  MappingTomorrow Jul 11 '12 at 7:18

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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