Python-GDAL-Tutorial: Kartenalgebra und Schreiben von Rasterdaten

Nehmen Sie die Berechnung von NDVI als Beispiel:

NDVI=(NIR-RED)/(NIR+RED)

wobei NIR Band 3 und RED Band 2 ist

Die Programmierpunkte lauten wie folgt:

1 Lesen Sie Band 3 in Array-Daten3 und lesen Sie Band 2 in Array-Daten2

2. Die Berechnungsformel lautet:

3. Wenn data3 und data2 beide 0 sind (0 steht beispielsweise für NODATA), tritt ein Division-durch-0-Fehler auf, der zum Absturz des Programms führt. Sie müssen mask mit „select“ verwenden, um den 0-Wert zu entfernen.

Der Code lautet wie folgt, es gibt nur 4 Zeilen

data2 = band2.ReadAsArray(0, 0, cols,rows). astype(Numeric.Float16)

data3 = band3.ReadAsArray(0, 0, cols,rows).astype(Numeric.Float16)

mask = Numeric.greater(data3 + data2, 0 )

ndvi = Numeric.choose(mask, (-99, (data3 - data2) / (data3 + data2)))

Neuen Rasterdatensatz erstellen

Schreiben die soeben berechneten Daten Im neuen Rasterdatensatz

kopieren Sie zunächst den Datentreiber:

driver = inDataset.GetDriver()

und erstellen Sie dann einen neuen Datensatz

Create(, , , [], [])

Der Standardwert von Bändern ist 1 und der Standardtyp von GDALDataType ist GDT_Byte, zum Beispiel

outDataset = drivers.Create(filename, cols, rows, 1, GDT_Float32)

Während der Ausführung dieser Anweisung wurde der Speicherplatz zugewiesen Festplatte

Vor dem Schreiben müssen Sie auch das Bandobjekt einführen

outBand = outDataset.GetRasterBand(1)

Das Bandobjekt unterstützt das direkte Schreiben in die Matrix und Die beiden Parameter sind der X-Richtungs-Offset und der Y-Richtungs-Offset

outBand.WriteArray(ndvi, 0, 0)

Das folgende Beispiel fasst die blockweise Schreibmethode dieses und letztes zusammen Zeit

xBlockSize = 64

yBlockSize = 64

für i in range(0, rows, yBlockSize):

if i + yBlockSize < :

numRows = yBlockSize

else:

numRows = rowsnumRows = rows –– ii

für j in range(0, cols, xBlockSize):

Wenn j + xBlockSize < cols:

numCols = xBlockSize

else:

numCols = cols – j

data = band.ReadAs Array( j, i, numCols, numRows)

# hier Berechnungen durchführen, um ein outData-Array zu erstellen

outBand.WriteArray(outData, j, i)

band Objekt kann NoData-Wert festlegen

outBand.SetNoDataValue(-99)

Sie können auch den NoData-Wert lesen

ND = outBand.GetNoDataValue()

Berechnen Sie die Statistiken des Bandes

Verwenden Sie zuerst FlushCache(), um Cache-Daten auf die Festplatte zu schreiben

Verwenden Sie dann GetStatistics(< approx_ok>, ), um Statistiken zu berechnen. Wenn approx_ok=1, basiert die Berechnung auf der Pyramide, wenn force=0, werden die Statistiken nicht berechnet, wenn das gesamte Diagramm erneut gelesen werden muss.

outBand.FlushCache()

outBand.GetStatistics(0, 1)

Legen Sie den geografischen Referenzpunkt des neuen Bildes fest

Wenn das neue Bild unterscheidet sich von anderen Wenn die geografischen Referenzinformationen einer Karte vollständig konsistent sind, ist es sehr einfach

geoTransform = inDataset.GetGeoTransform()

outDataset.SetGeoTransform(geoTransform)

proj = inDataset. GetProjection()

outDataset.SetProjection(proj)

Pyramiden erstellen

Pyramiden im Imagine-Stil festlegen

gdal.SetConfigOption('HFA_USE_RRD ', ' JA')

Pyramiden erzwingen

outDataset.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])

Zusammenfügen von Bildern

1. Für jedes Bild: Lesen Sie die Anzahl der Zeilen und Spalten, den Ursprung (minX, maxY), die Pixellänge und die Pixelbreite ab und berechnen Sie den Koordinatenbereich

maxX1 = minX1 + (cols1 * pixelWidth)

minY1 = maxY1 + (rows1 * pixelHeight)

2 Berechnen Sie den Koordinatenbereich des Ausgabebildes:

minX = min(minX1, minX2, … ) maxX = max( maxX1, maxX2, …)

minY = min(minY1, minY2, …) maxY = max(maxY1, maxY2, …)

3. Berechnen Sie die Anzahl der Zeilen und Spalten des Ausgabebildes:

cols = int((maxX – minX) / pixelWidth)

rows = int((maxY – minY) / abs(pixelHeight)

4. Erstellen und initialisieren Sie das Ausgabebild

5. Berechnen Sie den Offsetwert

xOffset1 = int((minX1 - minX) / pixelWidth)

yOffset1 = int((maxY1 - maxY) / pixelHeight)

Lesen Sie die Daten und schreiben Sie sie entsprechend dem oben berechneten Offset

6. Für das Ausgabebild: Statistiken berechnen und Geotransformation festlegen : [minX, pixelWidth, 0, maxY, 0, pixelHeight], Projektion festlegen, Pyramiden erstellen

Das Obige ist der Inhalt des Python-Gdal-Tutorials: Kartenalgebra und Rasterdatenschreiben. Weitere verwandte Inhalte finden Sie auf der chinesischen PHP-Website (www.php.cn)!