Tworzenie obrazu wielospektralnego od zera

10

Chcę zrobić obraz wielospektralny z cero, aby wykonać na nim kilka testów. Coś naprawdę prostego, jak 5 całkowicie jednolitych pasm z szumem soli i pieprzu lub kwadrat o różnych wartościach na środku. Najwyraźniej byłby to po prostu zestaw macierzy, wielowymiarowa tablica, którą generowanie jest dość proste. Chcę to osiągnąć za pomocą pytona i gdala, ale gdal jest dość hermetyczny, nie rozumiem tego wcale. Idealnie chciałbym utworzyć plik geotiff. Czy ktoś mógłby mi w tym pomóc? jakieś wskazówki lub tutorial gdal, który jest bardzo delikatny? Dziękuję wam wszystkim.

python raster gdal imagery JEquihua
źródło

Odpowiedzi:

15

Chcesz metodę gdal.band.WriteArray. W samouczku GDAL API znajduje się przykład (odtworzony poniżej):

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

Aby wygenerować losowe dane, spójrz na moduł numpy.random .

Oto bardziej kompletny działający przykład:

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
użytkownik2856
źródło
Wielkie dzięki, gdzie mogę przeczytać, co te rzeczy robią? SetUTM (ok, wiem, co to robi) SetWellKnown GeogCS, wyświetlanie, ustawianie transformacji geograficznej itp. ... ale wygląda dokładnie tak, jak potrzebuję. Wielkie dzięki!
JEquihua,
Więcej informacji na temat części kodu do georeferencji można znaleźć w samouczku Projekcje - gdal.org/ogr/osr_tutorial.html
user2856
2

Wiem, że nie o to prosiłeś, ale jeśli wszystko, czego potrzebujesz, to przykładowe dane wielospektralne lub hiperspektralne - te dane testowe dla projektu Opticks mogą działać. Alternatywnie możesz uzyskać dane LANDSAT bezpośrednio z Earth Explorer .

Ta strona zawiera przykładowy kod do konwersji dwuwymiarowej tablicy numpy na jednopasmowy geoTIFF, a wielopasmowy geoTIFF na 3D tablicę numpy.

EDYTOWAĆ:

Dalsze badania znajdują stronę przykładowego kodu z „brakującym przykładem”, tablicą numeryczną 3D -> wielopasmowym geoTIFF.

Mapowanie jutro
źródło
Nie, naprawdę muszę stworzyć swój własny wizerunek. Strona jest interesująca, dziękuję, tak naprawdę potrzebowałbym brakującego przykładu, jak zapisać tablicę 3d numpy jako wielopasmowy geoTIFF. Ale wielkie dzięki!
JEquihua,
Edytowano z dodatkowymi informacjami
MappingTomorrow