Criando uma imagem multiespectral do zero

Quero fazer uma imagem multiespectral de cero para fazer alguns testes nela. Algo realmente simples, como 5 bandas completamente uniformes com ruído de sal e pimenta ou um quadrado de valores diferentes no centro. Claramente, isso seria apenas uma pilha de matrizes, uma matriz multidimensional, que é bastante simples de gerar. Eu quero conseguir isso usando python e gdal, mas gdal está sendo bastante hermético, não entendo nada disso. Idealmente, eu gostaria de criar um arquivo geotiff. Alguém poderia me ajudar com isso? alguns ponteiros ou algum tutorial gdal que é muito gentil? Obrigado a todos.

Você deseja o método gdal.band.WriteArray. Há um exemplo no tutorial da API GDAL (reproduzido abaixo):

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

Para gerar os dados aleatórios, consulte o módulo numpy.random .

Aqui está um exemplo de trabalho mais completo:

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

Sei que não é o que você pediu, mas se tudo o que você deseja são dados de amostra multiespectrais ou hiperespectrais - esses dados de teste para o projeto Opticks podem funcionar. Como alternativa, você pode obter dados do LANDSAT diretamente do Earth Explorer .

Este site possui um exemplo de código para converter uma matriz numpy 2D em um geoTIFF de banda única e um geoTIFF multibanda em uma matriz numpy 3D.

EDITAR:

Pesquisas adicionais encontram uma página de código de exemplo com o 'exemplo ausente', matriz numpy 3D -> geoTIFF multi-banda.