Python-GDAL-Tutorial: Lesen von Rasterdaten mit GDAL

GDAL unterstützt nativ mehr als 100 Rasterdatentypen und deckt alle gängigen GIS- und RS-Datenformate ab, einschließlich

ArcInfo-Gitter, ArcSDE-Raster, Imagine, Idrisi, ENVI, GRASS, GeoTIFF

HDF4, HDF5

USGS DOQ, USGS DEM

ECW, MrSID

TIFF, JPEG, JPEG2000, PNG, GIF, BMP

Vollständige Unterstützung für Liste finden Sie unter http://www.gdal.org/formats_list.html

GDAL-Unterstützungsbibliothek importieren

Alte Version (vor 1.5): gdal, gdalconst importieren

Neue Version (1.6 und höher): aus Osgeo Import gdal, gdalconst

Es ist am besten, sowohl gdal als auch gdalconst zu importieren. Den Konstanten in gdalconst wird ein Präfix vorangestellt, um Konflikte mit anderen Modulen zu minimieren. So können Sie gdalconst direkt wie folgt importieren: from osgeo.gdalconst import *

Der GDAL-Datentreiber muss, ähnlich wie der OGR-Datentreiber, zuerst einen bestimmten Datentreibertyp erstellen und dann den entsprechenden Rasterdatensatz erstellen.

Registrieren Sie alle Datentreiber auf einmal, können aber nur lesen und nicht schreiben: gdal.AllRegister()

Registrieren Sie einen bestimmten Datentreibertyp separat, damit er gelesen und geschrieben werden kann. ja Erstellen Sie einen neuen Datensatz:

driver = gdal.GetDriverByName('HFA')

driver.Register()

Öffnen Sie einen vorhandenen Rasterdatensatz:

fn = 'aster.img'

ds = gdal.Open(fn, GA_ReadOnly)

wenn ds None ist:

print 'Konnte nicht geöffnet werden ' + fn

sys.exit(1)

Lesen Sie die Anzahl der Pixel in x-Richtung, die Anzahl der Pixel in y-Richtung und die Anzahl der Bänder des Rasterdatensatzes

cols = ds.RasterXSize

rows = ds.RasterYSize

bands = ds.RasterCount

Beachten Sie, dass danach keine Klammern stehen, da es sich um Eigenschaften (Eigenschaften) handelt. keine Methoden (Methoden)

Referenzinformationen zu geografischen Koordinaten (Georeferenzinformationen) lesen

GeoTransform ist eine Liste, die die geografischen Koordinateninformationen des Rasterdatensatzes speichert

adfGeoTransform[0] /* oben links x x-Koordinate der oberen linken Ecke */

adfGeoTransform[1] /* w--e Pixelauflösung Pixelauflösung in Ost-West-Richtung*/

adfGeoTransform[2 ] /* Rotation, 0, wenn das Bild „Norden oben“ ist. Wenn der Norden nach oben zeigt, der Rotationswinkel der Karte*/

adfGeoTransform[3] /* oben links y y-Koordinate der oberen linken Ecke* /

adfGeoTransform[4] /* Rotation, 0, wenn das Bild „nach Norden oben“ zeigt. Wenn der Norden nach oben zeigt, der Rotationswinkel der Karte*/

adfGeoTransform[5] /* n-s Pixelauflösung Pixelauflösung in Nord-Süd-Richtung*/

Beachten Sie die Rasterdaten. Die Koordinaten des Satzes basieren im Allgemeinen auf der oberen linken Ecke.

Das folgende Beispiel besteht darin, die Geotransformation aus einem Rasterdatensatz als Liste zu übernehmen und dann die darin enthaltenen Daten zu lesen

geotransform = ds.GetGeoTransform()

originX = geotransform [0]

originY = geotransform[3] originY = geotransform[3]

pixelWidth = geotransform[1]

pixelHeight = geotransform[5]

Berechnen Sie die relative Position (Pixelversatz) des Pixels, das einer bestimmten Koordinate entspricht, dh die relative Position der Koordinate und des Pixels in der oberen linken Ecke, berechnet basierend auf der Anzahl der Pixel. Die Berechnungsformel lautet wie folgt folgt:

xOffset = int((x – originX) / pixelWidth)

yOffset = int((y – originY) / pixelHeight)

Lesen des Werts eines bestimmten Pixel erfordert zwei Schritte

Zuerst lesen Holen Sie sich ein Band: GetRasterBand(), sein Parameter ist die Indexnummer des Bandes

und verwenden Sie dann ReadAsArray(, < ;yoff>, ,

band = ds.GetRasterBand(1)

data = band.ReadAsArray(xOffset, yOffset, 1, 1)

Wenn Sie ein Ganzes lesen möchten Band für Bild, dann setzen Sie den Offset auf 0 und die Größe auf die Größe des gesamten Bildes, zum Beispiel:

data = band.ReadAsArray(0, 0, cols, rows)

Bitte beachten Sie jedoch, dass Sie data[yoff, xoff] verwenden müssen, um den Wert eines bestimmten Pixels aus Daten zu lesen. Achten Sie darauf, es nicht rückwärts zu machen. Eine Matrix in der Mathematik ist [row,col], aber hier ist es genau das Gegenteil! Hier entspricht die Zeile der y-Achse und die Spalte der x-Achse.

Achten Sie darauf, bei Bedarf Speicher freizugeben, z. B. band = None oder dataset = None. Vor allem, wenn das Bild sehr groß ist

Wie liest man Rasterdaten effizienter? Offensichtlich ist das einzelne Lesen sehr ineffizient. Es ist keine gute Idee, den gesamten Rasterdatensatz in ein zweidimensionales Array zu stopfen, da dies immer noch viel Speicher beansprucht. Ein besserer Ansatz besteht darin, auf Daten in Blöcken zuzugreifen und nur den benötigten Block im Speicher abzulegen. Der Beispielcode dieser Woche verfügt über ein Utils-Modul, das die Blockgröße lesen kann.

Zum Beispiel:

Utils importieren

blockSize = utils.GetBlockSize(band)

xBlockSize = blockSize[0]

yBlockSize = blockSize[1]

Kacheln (gekachelt), das heißt, Rasterdaten werden in Blöcken gespeichert. Einige Formate, wie z. B. GeoTiff, sind nicht gekachelt und eine Zeile ist ein Block. Das Erdas-Imagine-Format ist mit 64x64 Pixeln gekachelt.

Wenn eine Zeile ein Block ist, ist das zeilenweise Lesen ressourcenschonender.

Wenn es sich um eine gekachelte Datenstruktur handelt, ist es die effizienteste Methode, den Parameterwert von ReadAsArray() so festzulegen, dass jeweils nur ein Block gelesen werden kann. Zum Beispiel:

rows = 13, cols = 11, xBSize = 5, yBSize = 5

for i in range(0, rows, yBSize):

if i + yBSize < rows:

numRows = yBSize

else:

numRows = rows – i

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

if j + xBSize < cols:

numCols = xBSize

else:

numCols = colsnumCol s = cols – j

data = band.ReadAsArray(j, i, numCols, numRows)

Dieser Code ist universell und kann von Zeit zu Zeit verwendet werden.

Im Folgenden werden einige Fähigkeiten zur zweidimensionalen Array-Verarbeitung vorgestellt.

Hier werden zwei Bibliotheken verwendet: Numerisch und Numpy. Numerisch ist älter und wird von FWTools verwendet. Wenn Sie es selbst installieren und konfigurieren, sollten Sie das leistungsfähigere Numpy verwenden.

Datentypkonvertierung:

data = band.ReadAsArray(j, i, nCols, nRows)

data = data.astype(Numeric.Float) # Numeric

data = data.astype(numpy.float) # numpy

Oder schreiben Sie einfach einen einfachen Satz

data = band.ReadAsArray(j, i, nCols, nRows).astype (Numeric.Float)

Maskenmaske

Dies ist die Funktion von Numeric- und Numpy-Bibliotheken. Sie gibt ein Array und eine Bedingung ein und gibt ein binäres Array aus. Zum Beispiel

mask = Numeric.greater(data, 0)mask = Numeric.greater(data, 0)

>>> , 6, 0, 2])

>>> print a

[0 4 6 0 2]

>>> mask = Numerisch . great(a, 0)

>>> print mask

[0 1 1 0 1]

Array-Summierung

> >> a = Numeric.array([0, 4, 6, 0, 2])

>>> print a

0 4 6 0 2]

>>> print Numeric.sum(a)

12

Wenn es sich um ein zweidimensionales Array handelt, dann wird sum ein eindimensionales Array zurückgeben

>>> b = Numeric.array([a, [5, 10, 0, 3, 0]])

>> > print b

[[ 0 4 6 0 2]

[ 5 10 0 3 0]]

>>> print Numeric.sum(b )> >> print Numeric.sum(b)

[ 5 14 6 3 2]

Die Summe zweidimensionaler Arrays ist also so

>> ;> print Numeric.sum(Numeric.sum(b))

30

Hier ist ein kleiner Trick, um die Anzahl der Pixel größer als 0 zu zählen. Sie können verwenden Maske und Summe zusammen. Funktion

>>> print a

[0 4 6 0 2]

>>> (a, 0)

>>> print mask

[0 1 1 0 1]

>>> print Numeric.sum(mask )

3

Das Obige ist das Python-Gdal-Tutorial: Verwenden von Gdal zum Lesen von Rasterdaten. Weitere verwandte Inhalte finden Sie auf der chinesischen PHP-Website (www.php.cn)!