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?
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 nrows=1024 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 upper left corner of the image #in same units as spatial reference #xmin=147.2 #ymax=-34.54 #Cellsize in same units as spatial reference #cellsize=0.01 #dst_ds.SetGeoTransform( [ xmin, cellsize, 0, ymax, 0, -cellsize ] ) #srs = osr.SpatialReference() #srs.SetWellKnownGeogCS("WGS84") #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
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.