GDAL und Python: Wie erhält man Koordinaten für alle Zellen mit einem bestimmten Wert?

Ich habe ein Arc / Info-Binärgitter - speziell ein ArcGIS-Flow-Akkumulations-Raster - und möchte alle Zellen mit einem bestimmten Wert (oder einem Wertebereich) identifizieren. Letztendlich hätte ich gerne ein Shapefile mit Punkten, die diese Zellen darstellen.

Ich kann QGIS verwenden, um die hdr.adf zu öffnen und dieses Ergebnis zu erhalten. Der Workflow ist:

• QGIS> Menü Raster> Rasterrechner (alle Punkte mit Zielwert markieren)
• QGIS> Raster-Menü> Polygonisieren
• QGIS> Menü Vektor> Untermenü Geometrie> Polygon-Schwerpunkte
• Bearbeiten Sie die Centroids, um die unerwünschten Polycentroids zu löschen (those = 0).

Dieser Ansatz "erledigt den Job", spricht mich aber nicht an, da er 2 Dateien erstellt, die ich löschen muss, und dann unerwünschte Datensätze aus dem Shapefile von Zentroiden (dh jene = 0) entfernen muss.

Eine vorhandene Frage nähert sich diesem Thema, ist jedoch auf ArcGIS / ArcPy zugeschnitten, und ich möchte im FOSS-Bereich bleiben.

Verfügt jemand über ein vorhandenes GDAL / Python-Rezept / -Skript, das die Zellwerte eines Rasters abfragt, und wenn ein Zielwert --- oder ein Wert in einem Zielbereich --- gefunden wird, wird einem Shapefile ein Datensatz hinzugefügt? Dies würde nicht nur die Interaktion mit der Benutzeroberfläche vermeiden, sondern in einem einzigen Durchgang ein sauberes Ergebnis erzielen.

Ich habe einen Versuch unternommen, indem ich gegen eine der Präsentationen von Chris Garrard gearbeitet habe , aber die Rasterarbeit befindet sich nicht in meinem Steuerhaus und ich möchte die Frage nicht mit meinem schwachen Code überfrachten.

Wenn jemand den genauen Datensatz zum Spielen haben möchte, füge ich ihn hier als .zip ein .

[Notizen bearbeiten] Lassen Sie dies für die Nachwelt hinter sich. Siehe Kommentaraustausch mit om_henners. Grundsätzlich wurden die x / y-Werte (Zeile / Spalte) gespiegelt. Die ursprüngliche Antwort hatte diese Zeile:

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

umgekehrt, so:

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

Als ich zum ersten Mal auf das im Screenshot dargestellte Problem stieß, fragte ich mich, ob ich die Geometrie falsch implementiert hatte, und experimentierte, indem ich die x / y-Koordinatenwerte in dieser Zeile spiegelte:

point.SetPoint(0, x, y)

..wie..

point.SetPoint(0, y, x)

Das hat aber nicht funktioniert. Und ich dachte nicht daran, die Werte im Numpy-Ausdruck von om_henners zu spiegeln, weil ich fälschlicherweise glaubte, dass das Spiegeln in beiden Zeilen gleichwertig ist. Ich denke , das eigentliche Problem bezieht sich auf die x_sizeund y_sizeWerte, jeweils 30und -30, die angewendet werden , wenn die Zeilen- und Spaltenindizes zu berechnen Punktkoordinaten für die Zellen verwendet werden.

[Original Edit]

@om_henners, ich versuche Ihre Lösung zusammen mit ein paar Rezepten zum Erstellen von Punkt-Shapefiles mit ogr ( invisibleroads.com , Chris Garrard ), aber ich habe ein Problem, bei dem die Punkte so erscheinen, als würden sie sich über eine Linie hinweg widerspiegeln durch 315/135-Grad.

Hellblaue Punkte : Erstellt mit meinem QGIS-Ansatz oben

Lila Punkte : Erstellt mit dem unten stehenden GDAL / OGR-PY-Code

[Gelöst]

Dieser Python-Code implementiert die von @om_henners vorgeschlagene Komplettlösung. Ich habe es getestet und es funktioniert. Danke, Mann!

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!"

In GDAL können Sie das Raster als Numpy-Array importieren.

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)

Dann ist es mit numpy ganz einfach, die Indizes eines Arrays abzurufen, die mit einem boolan-Ausdruck übereinstimmen:

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

Aus der Raster-Geotransformation können Informationen wie die oberen linken x- und y-Koordinaten und die Zellengrößen abgerufen werden.

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

Die obere linke Zelle entspricht a[0, 0]. Die Y-Größe ist immer negativ. Mit den X- und Y-Indizes können Sie die Koordinaten jeder Zelle basierend auf den Indizes berechnen.

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

Ab hier ist es ganz einfach, mit OGR ein Shapefile zu erstellen. Beispielcode finden Sie in dieser Frage zum Generieren eines neuen Datasets mit Punktinformationen.

Warum nicht die Sexante-Toolbox in QGIS verwenden? Es ähnelt dem Model Builder für ArcGIS. Auf diese Weise können Sie Vorgänge verketten und als einen Vorgang behandeln. Sie können das Löschen von Zwischendateien und das Entfernen unerwünschter Datensätze automatisieren, wenn ich mich nicht irre.

Es kann hilfreich sein, die Daten nach Postgis zu importieren (mit Rasterunterstützung) und die dortigen Funktionen zu verwenden. Dieses Tutorial enthält möglicherweise Elemente, die Sie benötigen.