Criando uma imagem multiespectral do zero

10

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.

JEquihua
fonte

Respostas:

15

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
user2856
fonte
Muito obrigado, onde podemos ler o que essas coisas fazem? SetUTM (ok, eu sei o que isso faz) SetWellKnown GeogCS, se projeção, set geo transform, etc ... mas parece exatamente o que eu preciso. Muito obrigado!
JEquihua
Para mais informações sobre as partes de georreferenciamento do código, consulte o Tutorial de projeções - gdal.org/ogr/osr_tutorial.html
user2856
2

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.

MappingAmanhã
fonte
Não, eu realmente preciso criar minha própria imagem. A página é interessante, obrigado, o que eu realmente precisaria é o exemplo que faltava, como salvar uma matriz numpy 3d como um geoTIFF de várias bandas. Mas muito obrigado!
JEquihua
Editado com mais informações
MappingTomorrow