Membuat gambar multispektral dari awal

10

Saya ingin membuat gambar multispektral dari cero untuk melakukan beberapa tes di atasnya. Sesuatu yang sangat sederhana seperti 5 band yang sepenuhnya seragam dengan bunyi garam dan merica di atasnya atau kuadrat dari nilai yang berbeda di tengah. Jelas ini hanya akan menjadi tumpukan matriks, array multidimensi, yang cukup mudah untuk dihasilkan. Saya ingin mencapai ini menggunakan python dan gdal tetapi gdal sedang cukup hermetis, saya tidak mengerti sama sekali. Saya idealnya ingin membuat file geotiff. Adakah yang bisa membantu saya dalam hal ini? beberapa petunjuk atau beberapa tutorial gdal yang sangat lembut? Terima kasih semua.

python raster gdal imagery JEquihua
sumber

Jawaban:

15

Anda menginginkan metode gdal.band.WriteArray. Ada contoh dalam tutorial API GDAL (direproduksi di bawah):

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

Untuk menghasilkan data acak, lihat modul numpy.random .

Berikut ini contoh kerja yang lebih lengkap:

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
pengguna2856
sumber
Terima kasih banyak, di mana bisa membaca apa yang dilakukan hal-hal ini? SetUTM (ok saya tahu apa yang dilakukannya) SetWellKnown GeogCS, proyeksi se, atur geo transform, dll ... tetapi sepertinya persis seperti yang saya butuhkan. Terima kasih banyak!
JEquihua
Untuk info lebih lanjut tentang bagian-bagian georeferensi kode, lihat Tutorial Proyeksi - gdal.org/ogr/osr_tutorial.html
user2856
2

Saya tahu itu bukan yang Anda minta, tetapi jika yang Anda inginkan hanyalah data sampel multispektral atau hiperspektral - data uji ini untuk proyek Opticks mungkin berfungsi. Sebagai alternatif, Anda bisa mendapatkan data LANDSAT langsung dari Earth Explorer .

Situs ini memiliki kode contoh untuk mengubah larik numpy 2D ke geoTIFF pita tunggal, dan geoTIFF multi-band menjadi larik numpy 3D.

EDIT:

Penelitian lebih lanjut menemukan halaman kode contoh dengan 'contoh hilang', array numpy 3D -> multi-band geoTIFF.

Memetakan Besok
sumber
Tidak, saya benar-benar perlu membuat gambar saya sendiri. Halaman ini menarik, terima kasih, yang benar-benar saya butuhkan adalah contoh yang hilang, bagaimana cara menyimpan array numpy 3d sebagai multi-band geoTIFF. Tapi terima kasih banyak!
JEquihua
Diedit dengan lebih banyak info
MappingTomorrow