GDAL i Python: Jak uzyskać współrzędne dla wszystkich komórek mających określoną wartość?

Mam siatkę binarną Arc / Info --- konkretnie raster akumulacji przepływu ArcGIS --- i chciałbym zidentyfikować wszystkie komórki mające określoną wartość (lub zakres wartości). Ostatecznie chciałbym uzyskać plik kształtu reprezentujący te komórki.

Mogę użyć QGIS, aby otworzyć plik hdr.adf i uzyskać ten wynik, przepływ pracy jest następujący:

• QGIS> Menu rastrowe> Kalkulator rastrowy (zaznacz wszystkie punkty wartością docelową)
• QGIS> Menu rastrowe> Poligonizuj
• QGIS> Menu wektorowe> podmenu Geometria> Centroidy wielokątów
• Edytuj centroidy, aby usunąć niechciane centroidy (te = 0)

To podejście „wykonuje pracę”, ale nie podoba mi się, ponieważ tworzy 2 pliki, które muszę usunąć, a następnie muszę usunąć niechciane zapisy z pliku kształtu centroidów (tj. Te = 0).

Istniejące pytanie podchodzi do tego tematu, ale to jest dostosowane do ArcGIS / ArcPy i chciałbym zatrzymać się w przestrzeni FOSS.

Czy ktoś ma istniejącą recepturę / skrypt GDAL / Python, który przesłuchuje wartości komórek rastra, a gdy znaleziona zostanie wartość docelowa --- lub wartość z zakresu docelowego ---, rekord zostanie dodany do pliku kształtu? Pozwoli to nie tylko uniknąć interakcji z interfejsem użytkownika, ale zapewni czysty wynik w jednym przejściu.

Strzeliłem do niego, pracując przeciwko jednej z prezentacji Chrisa Garrarda , ale praca nad rastrami nie jest w mojej sterówce i nie chcę zaśmiecać tego pytania słabym kodem.

Jeśli ktoś chce grać z dokładnym zestawem danych, umieszczam go tutaj jako .zip .

[Edytuj notatki] Pozostawiając to dla potomności. Zobacz wymiany komentarzy z om_henners. Zasadniczo wartości x / y (wiersz / kolumna) zostały odwrócone. Oryginalna odpowiedź miała następującą linię:

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

odwrócony, jak poniżej:

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

Kiedy po raz pierwszy zetknąłem się z problemem przedstawionym na zrzucie ekranu, zastanawiałem się, czy nieprawidłowo zaimplementowałem geometrię, i eksperymentowałem, odwracając wartości współrzędnych x / y w tym wierszu:

point.SetPoint(0, x, y)

..tak jak..

point.SetPoint(0, y, x)

Jednak to nie zadziałało. I nie pomyślałem, żeby spróbować przerzucić wartości w wyrażeniu Numpy om_henners, niesłusznie wierząc, że odwrócenie ich w dowolnej linii jest równoważne. Myślę prawdziwy problem dotyczy x_sizeoraz y_sizewartości, odpowiednio, 30i -30, które są stosowane w przypadku gdy wierszy i kolumn Wskaźniki są wykorzystywane do obliczenia współrzędnych punktów dla komórek.

[Edycja oryginalna]

@om_henners, wypróbowuję twoje rozwiązanie w porozumieniu z kilkoma przepisami na tworzenie plików kształtów punktów za pomocą ogr ( invisibleroads.com , Chris Garrard ), ale mam problem, w którym punkty pojawiają się, jakby odbijały się w poprzek linii do 315/135 stopni.

Jasnoniebieskie punkty : utworzone przez moje podejście QGIS powyżej

Punkty fioletowe : utworzone poniżej za pomocą kodu py GDAL / OGR

[Rozwiązany]

Ten kod Python implementuje kompletne rozwiązanie zaproponowane przez @om_henners. Przetestowałem to i działa. Dzięki!

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

W GDAL możesz zaimportować raster jako tablicę 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)

Następnie za pomocą numpy proste jest uzyskanie indeksów tablicy pasujących do wyrażenia boolan:

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

Z geotransformy rastrowej możemy uzyskać takie informacje, jak lewe górne współrzędne xiy oraz rozmiary komórek.

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

Lewa górna komórka odpowiada a[0, 0]. Rozmiar Y zawsze będzie ujemny, więc używając indeksów xiy możesz obliczyć współrzędne każdej komórki na podstawie indeksów.

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

Odtąd wystarczy prosta kwestia, aby utworzyć plik kształtu za pomocą OGR. Aby uzyskać przykładowy kod, zobacz to pytanie, jak wygenerować nowy zestaw danych z informacjami o punkcie.

Dlaczego nie skorzystać z przybornika Sexante w QGIS? To trochę jak Konstruktor modeli dla ArcGIS. W ten sposób możesz połączyć operacje i traktować je jako jedną operację. Możesz zautomatyzować usuwanie plików pośrednich i usuwanie niechcianych rekordów, jeśli się nie mylę.

Pomocne może być zaimportowanie danych do Postgis (z obsługą rastra) i korzystanie z dostępnych tam funkcji. Ten samouczek może zawierać elementy, których potrzebujesz.