tutoriel python gdal : algèbre cartographique et écriture de données raster

Prenons comme exemple le calcul du NDVI :

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

où NIR est la bande 3 et RED est la bande 2

Programmation Les points clés sont les suivants :

1. Lisez la bande 3 dans les données du tableau3 et lisez la bande 2 dans les données du tableau2

2. La formule de calcul est :

3. Lorsque data3 et data2 Lorsque les deux valent 0 (par exemple, 0 représente NODATA), une erreur de division par 0 se produira, provoquant le crash du programme. Vous devez utiliser un masque avec choisir pour supprimer la valeur 0

Le code est le suivant, il n'y a que 4 lignes

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)))

Créer un nouvel ensemble de données raster

Écrire les données vient d'être calculé dans le nouveau Dans l'ensemble de données raster

copiez d'abord le pilote de données :

driver = inDataset.GetDriver()

puis créez un nouvel ensemble de données

Créer (, , , [], [])

La valeur par défaut des bandes est 1, et la valeur par défaut le type de GDALDataType est GDT_Byte , par exemple

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

Lors de l'exécution de cette instruction, l'espace de stockage a été alloué à le disque dur

Avant d'écrire, vous devez d'abord introduire l'objet bande

outBand = outDataset.GetRasterBand(1)

L'objet bande prend en charge l'écriture directe sur la matrice, et les deux paramètres sont le décalage dans la direction x et le décalage dans la direction y

outBand.WriteArray(ndvi, 0, 0)

L'exemple suivant résume cette méthode et la dernière méthode d'écriture bloc par bloc

xBlockSize = 64

yBlockSize = 64

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

si je yBlockSize < 🎜>

numRows = yBlockSize

else:

numRows = rowsnumRows = rows –– ii

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

si j xBlockSize < cols :

numCols = xBlockSize

else :

numCols = cols – j

data = band.ReadAsArray (j, i, numcols, numrows)

# Effectuez des calculs ici pour créer un tableau Outdata

Outband.writearray (OUTDATA, J, I)

La valeur Nodata

outBand.SetNoDataValue(-99)

Vous pouvez également lire la valeur NoData

ND = outBand.GetNoDataValue()

Calculer les statistiques du groupe

Utilisez d'abord FlushCache() pour écrire les données du cache sur le disque

, puis utilisez GetStatistics(, ) pour calculer les statistiques. Si approx_ok=1 alors le calcul est basé sur la pyramide, si force=0 alors les statistiques ne sont pas calculées lorsqu'il faut relire tout le graphique.

outBand.FlushCache()

outBand.GetStatistics(0, 1)

Définir le point de référence géographique de la nouvelle image

Si la nouvelle image est différent d'une autre Si les informations de référence géographique d'une carte sont tout à fait cohérentes, c'est très simple

geoTransform = inDataset.GetGeoTransform()

outDataset.SetGeoTransform(geoTransform)

proj = inDataset. GetProjection()

outDataset.SetProjection(proj)

Créer des pyramides

Définir les pyramides de style Imagine

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

Forcer la construction de pyramides

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

Assemblage d'images

1. Pour chaque image : lisez le nombre de lignes et de colonnes, l'origine (minX, maxY), la longueur des pixels, la largeur des pixels et calculez la plage de coordonnées

maxX1 = minX1 (cols1 * pixelWidth)

minY1 = maxY1 (rows1 * pixelHeight)

2 Calculez la plage de coordonnées de l'image de sortie :

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

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

3. Calculez le nombre de lignes et de colonnes. de l'image de sortie :

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

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

4 .         Créez et initialisez l'image de sortie

5. Pour chaque image à coller : Calculez la valeur de décalage

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

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

Lisez les données et écrivez-les en fonction du décalage calculé ci-dessus

6 Pour l'image de sortie : Calculez les statistiques et définissez la géotransformation : [ minX, pixelWidth, 0, maxY, 0, pixelHeight], définir la projection, créer des pyramides

Ce qui précède est le contenu du didacticiel python gdal : algèbre cartographique et écriture de données raster. Pour plus de contenu connexe, veuillez faire attention au site Web PHP chinois (www.php.cn) !