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!"
srs.ImportFromWkt(r.GetProjection())(zamiast konieczności tworzenia projekcji ze znanego ciągu proj).