GDAL et Python: Comment obtenir les coordonnées de toutes les cellules ayant une valeur spécifique?

J'ai une grille binaire Arc / Info --- spécifiquement, un raster d'accumulation de flux ArcGIS --- et je voudrais identifier toutes les cellules ayant une valeur spécifique (ou dans une plage de valeurs). En fin de compte, je voudrais un fichier de formes de points représentant ces cellules.

Je peux utiliser QGIS pour ouvrir le hdr.adf et obtenir ce résultat, le workflow est:

• QGIS> menu Raster> calculatrice raster (marquer tous les points avec la valeur cible)
• QGIS> Menu raster> Polygoniser
• QGIS> Menu vectoriel> Sous-menu Géométrie> Centroïdes de polygone
• Modifier les centroïdes pour supprimer les poly centroïdes indésirables (ceux = 0)

Cette approche "fait le travail", mais elle ne m'intéresse pas car elle crée 2 fichiers que je dois supprimer, puis je dois supprimer les enregistrements indésirables du fichier de formes des centroïdes (c'est-à-dire ceux = 0).

Une question existante aborde ce sujet, mais elle est adaptée à ArcGIS / ArcPy, et j'aimerais rester dans l'espace FOSS.

Quelqu'un a-t-il une recette / un script GDAL / Python existant qui interroge les valeurs des cellules d'un raster, et lorsqu'une valeur cible --- ou une valeur dans une plage cible --- est trouvée, un enregistrement est ajouté à un fichier de formes? Cela éviterait non seulement l'interaction avec l'interface utilisateur, mais créerait un résultat net en une seule passe.

J'ai pris une chance en travaillant contre l'une des présentations de Chris Garrard , mais le raster n'est pas dans ma timonerie et je ne veux pas encombrer la question avec mon code faible.

Si quelqu'un veut jouer avec l'ensemble de données exact, je le mets ici sous forme de .zip .

[Modifier les notes] Laissant cela derrière pour la postérité. Voir les échanges de commentaires avec om_henners. Fondamentalement, les valeurs x / y (ligne / colonne) ont été inversées. La réponse originale avait cette ligne:

(y_index, x_index) = np.nonzero(a == 1000)

inversé, comme ceci:

(x_index, y_index) = np.nonzero(a == 1000)

Lorsque j'ai rencontré le problème illustré dans la capture d'écran pour la première fois, je me suis demandé si j'avais mal mis en œuvre la géométrie et j'ai expérimenté en inversant les valeurs des coordonnées x / y dans cette ligne:

point.SetPoint(0, x, y)

..comme..

point.SetPoint(0, y, x)

Mais cela n'a pas fonctionné. Et je n'ai pas pensé à essayer de retourner les valeurs dans l'expression Numpy d'om_henners, croyant à tort que les retourner sur l'une ou l'autre ligne était équivalent. Je pense que le vrai problème concerne les valeurs x_sizeet y_size, respectivement, 30et -30, qui sont appliquées lorsque les indices de ligne et de colonne sont utilisés pour calculer les coordonnées des points pour les cellules.

[Édition originale]

@om_henners, j'essaie votre solution, de concert avec quelques recettes pour créer des fichiers de formes de points à l'aide d'ogr ( invisibleroads.com , Chris Garrard ), mais j'ai un problème où les points apparaissent comme s'ils se reflétaient sur une ligne passant à 315/135 degrés.

Points bleu clair : créés par mon approche QGIS , ci-dessus

Points violets : créés par le code py GDAL / OGR , ci-dessous

[Résolu]

Ce code Python implémente la solution complète proposée par @om_henners. Je l'ai testé et ça marche. Merci mec!

from osgeo import gdal
import numpy as np
import osgeo.ogr
import osgeo.osr

path = "D:/GIS/greeneCty/Greene_DEM/GreeneDEM30m/flowacc_gree/hdr.adf"
print "\nOpening: " + path + "\n"

r = gdal.Open(path)
band = r.GetRasterBand(1)

(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = r.GetGeoTransform()

a = band.ReadAsArray().astype(np.float)

# This evaluation makes x/y arrays for all cell values in a range.
# I knew how many points I should get for ==1000 and wanted to test it.
(y_index, x_index) = np.nonzero((a > 999) & (a < 1001))

# This evaluation makes x/y arrays for all cells having the fixed value, 1000.
#(y_index, x_index) = np.nonzero(a == 1000)

# DEBUG: take a look at the arrays..
#print repr((y_index, x_index))

# Init the shapefile stuff..
srs = osgeo.osr.SpatialReference()
#srs.ImportFromProj4('+proj=utm +zone=15 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')
srs.ImportFromWkt(r.GetProjection())

driver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
shapeData = driver.CreateDataSource('D:/GIS/01_tutorials/flow_acc/ogr_pts.shp')

layer = shapeData.CreateLayer('ogr_pts', srs, osgeo.ogr.wkbPoint)
layerDefinition = layer.GetLayerDefn()

# Iterate over the Numpy points..
i = 0
for x_coord in x_index:
    x = x_index[i] * x_size + upper_left_x + (x_size / 2) #add half the cell size
    y = y_index[i] * y_size + upper_left_y + (y_size / 2) #to centre the point

    # DEBUG: take a look at the coords..
    #print "Coords: " + str(x) + ", " + str(y)

    point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
    point.SetPoint(0, x, y)

    feature = osgeo.ogr.Feature(layerDefinition)
    feature.SetGeometry(point)
    feature.SetFID(i)

    layer.CreateFeature(feature)

    i += 1

shapeData.Destroy()

print "done! " + str(i) + " points found!"

Dans GDAL, vous pouvez importer le raster en tant que tableau numpy.

from osgeo import gdal
import numpy as np

r = gdal.Open("path/to/raster")
band = r.GetRasterBand(1) #bands start at one
a = band.ReadAsArray().astype(np.float)

Ensuite, en utilisant numpy, il est simple d'obtenir les index d'un tableau correspondant à une expression booléenne:

(y_index, x_index) = np.nonzero(a > threshold)
#To demonstate this compare a.shape to band.XSize and band.YSize

À partir de la géotransformation raster, nous pouvons obtenir des informations telles que les coordonnées x et y en haut à gauche et la taille des cellules.

(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = r.GetGeoTransform()

La cellule supérieure gauche correspond à a[0, 0]. La taille Y sera toujours négative, donc en utilisant les indices x et y, vous pouvez calculer les coordonnées de chaque cellule en fonction des indices.

x_coords = x_index * x_size + upper_left_x + (x_size / 2) #add half the cell size
y_coords = y_index * y_size + upper_left_y + (y_size / 2) #to centre the point

À partir d'ici, il suffit de créer un fichier de formes à l'aide d'OGR. Pour un exemple de code, consultez cette question pour savoir comment générer un nouvel ensemble de données avec des informations de point.

Pourquoi ne pas utiliser la boîte à outils Sexante dans QGIS? C'est un peu comme le Model Builder pour ArcGIS. De cette façon, vous pouvez chaîner des opérations et les traiter comme une seule opération. Vous pouvez automatiser la suppression des fichiers intermédiaires et la suppression des enregistrements indésirables si je ne me trompe pas.

Il peut être utile d'importer les données dans postgis (avec prise en charge raster) et d'utiliser les fonctions qui s'y trouvent. Ce didacticiel peut contenir des éléments dont vous avez besoin.