tutoriel python gdal : lecture de données raster avec gdal

GDAL prend en charge nativement plus de 100 types de données raster, couvrant tous les formats de données SIG et RS courants, y compris

grilles ArcInfo, raster ArcSDE, Imagine, Idrisi, ENVI, GRASS, GeoTIFF

HDF4, HDF5

USGS DOQ, USGS DEM

ECW, MrSID

TIFF, JPEG, JPEG2000, PNG, GIF, BMP

Support complet pour le liste, veuillez vous référer à http://www.gdal.org/formats_list.html

Importer la bibliothèque de support GDAL

Ancienne version (avant 1.5) : importer gdal, gdalconst

Nouvelle version (1.6 et versions ultérieures) : depuis osgeo import gdal, gdalconst

Il est préférable d'importer à la fois gdal et gdalconst. Les constantes de gdalconst sont préfixées pour minimiser les conflits avec d'autres modules. Vous pouvez donc importer directement gdalconst comme ceci : depuis osgeo.gdalconst import *

Le pilote de données GDAL, similaire au pilote de données OGR, doit d'abord créer un certain type de pilote de données, puis créer l'ensemble de données raster correspondant.

Enregistrez tous les pilotes de données à la fois, mais ne pouvez que lire et non écrire : gdal.AllRegister()

Enregistrez un certain type de pilote de données séparément, afin qu'il puisse être lu et écrit, oui Créer un nouveau jeu de données :

driver = gdal.GetDriverByName('HFA')

driver.Register()

Ouvrir un jeu de données raster existant :

fn = 'aster.img'

ds = gdal.Open(fn, GA_ReadOnly)

si ds est Aucun :

imprimer 'Impossible d'ouvrir ' fn

sys.exit(1)

Lire le nombre de pixels dans la direction x, le nombre de pixels dans la direction y et le nombre de canaux de l'ensemble de données raster

cols = ds.RasterXSize

rows = ds.RasterYSize

bands = ds.RasterCount

Notez qu'il n'y a pas de crochets après, car ce sont des propriétés (propriétés) non méthodes (méthodes)

Lire Obtenir des informations de référence de coordonnées géographiques (informations de géoréférence)

GeoTransform est une liste qui stocke les informations de coordonnées géographiques de l'ensemble de données raster

adfGeoTransform[0] /* coin supérieur gauche x coin supérieur gauche x coordonnée* /

adfGeoTransform[1] /* w--e résolution en pixels résolution en pixels dans la direction est-ouest*/

adfGeoTransform[2] /* rotation, 0 si l'image est "nord vers le haut" Si le nord est orienté vers le haut, l'angle de rotation de la carte*/

adfGeoTransform[3] /* en haut à gauche y y coordonnée du coin supérieur gauche*/

adfGeoTransform[4] /* rotation, 0 si l'image est "nord vers le haut" Si le nord est vers le haut, l'angle de rotation de la carte*/

adfGeoTransform[5] /* n-s résolution en pixels résolution en pixels dans la direction nord-sud*/

Notez le jeu de données raster. Les coordonnées sont généralement basées sur le coin supérieur gauche.

L'exemple suivant consiste à prendre le Geotransform d'un ensemble de données raster sous forme de liste, puis à lire les données qu'il contient

geotransform = ds.GetGeoTransform()

originX = géotransformation [0]

originY = géotransform[3] originY = géotransform[3]

pixelWidth = géotransform[1]

pixelHeight = géotransformation[5]

Calculez la position relative (décalage de pixels) du pixel correspondant à une certaine coordonnée, c'est-à-dire la position relative de la coordonnée et du pixel dans le coin supérieur gauche, calculée en fonction du nombre de pixels. La formule de calcul est la suivante. suit :

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

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

Lecture de la valeur d'un certain le pixel nécessite deux étapes

Lire d'abord Obtenir une bande : GetRasterBand(), son paramètre est le numéro d'index de la bande

puis utiliser ReadAsArray(, < ;yoff>, ,

band = ds.GetRasterBand(1)

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

Si vous souhaitez lire un texte entier bande à une image temporelle, puis définissez le décalage sur 0 et la taille sur la taille de l'image entière, par exemple :

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

Mais veuillez noter que pour lire la valeur d'un certain pixel à partir de données, vous devez utiliser data[yoff, xoff]. Attention à ne pas le faire à l'envers. Une matrice en mathématiques c’est [row,col], mais ici c’est tout le contraire ! Ici, row correspond à l'axe y et col correspond à l'axe x.

Faites attention à libérer de la mémoire lorsque cela est approprié, comme band = None ou dataset = None. Surtout lorsque l'image est très grande

Comment lire les données raster plus efficacement ? Évidemment, lire un par un est très inefficace. Ce n'est pas une bonne idée de regrouper l'intégralité du jeu de données raster dans un tableau à deux dimensions, car cela prend encore beaucoup de mémoire. Une meilleure approche consiste à accéder aux données par blocs et à mettre uniquement le bloc dont vous avez besoin en mémoire. L'exemple de code de cette semaine comporte un module utils capable de lire la taille du bloc.

Par exemple :

import utils

blockSize = utils.GetBlockSize(band)

xBlockSize = blockSize[0]

yBlockSize = taille de bloc[1]

Tile (carrelé), c'est-à-dire que les données raster sont stockées en blocs. Certains formats, comme GeoTiff, ne sont pas carrelés et chaque ligne est un bloc. Le format Erdas imagine est carrelé de 64x64 pixels.

Si une ligne est un bloc, alors la lecture par ligne permet d'économiser davantage de ressources.

S'il s'agit d'une structure de données en mosaïque, définir la valeur du paramètre ReadAsArray() afin qu'il ne puisse lire qu'un seul bloc à la fois est la méthode la plus efficace. Par exemple :

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

pour i dans la plage (0, lignes, yBSize) :

si je yBSize < lignes :

numRows = yBSize

else :

numRows = lignes – i

pour j dans la plage (0, cols, xBSize) :

si j xBSize < cols :

numCols = xBSize

sinon :

numCols = colsnumCols = col s – j

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

Ce code est universel et peut être utilisé de temps en temps.

Ce qui suit présente quelques compétences en matière de traitement de tableaux bidimensionnels

Deux bibliothèques sont utilisées ici, Numeric et numpy. Le numérique est plus ancien et FWTools l'utilise. Si vous l'installez et le configurez vous-même, vous devez utiliser le numpy plus puissant.

Conversion de type de données :

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

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

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

Ou écrivez simplement une phrase simple

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

Masque

C'est la fonction des bibliothèques numériques et numpy. Il entre un tableau et une condition et génère un tableau binaire. Par exemple

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

>>> , 6, 0, 2])

>>> imprimer un

[0 4 6 0 2]

>>> . Greater(a, 0)

>>> masque d'impression

[0 1 1 0 1]

Somme du tableau

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

>>> imprimer a>> imprimer un

[ 0 4 6 0 2]

>>> print Numeric.sum(a)

12

S'il s'agit d'un tableau à deux dimensions, alors la somme sera renvoie un tableau unidimensionnel

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

>> > imprimer b

[[ 0 4 6 0 2]

[ 5 10 0 3 0]]

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

[ 5 14 6 3 2]

Donc, la somme des tableaux bidimensionnels est comme ça

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

30

Voici une petite astuce Pour compter le nombre de pixels supérieur à 0, vous pouvez. utiliser le masque et la somme ensemble. Fonction

>>> imprimer un

[0 4 6 0 2]

>>> supérieur(a, 0)

>>> masque d'impression

[0 1 1 0 1]

>>> masque)

3

Ce qui précède est le tutoriel python gdal : utiliser gdal pour lire des données raster. Pour plus de contenu connexe, veuillez faire attention au site Web PHP chinois (www.php.cn). !