GDAL dan Python: Bagaimana cara mendapatkan koordinat untuk semua sel yang memiliki nilai tertentu?

Saya punya Arc / Info Binary Grid --- khususnya, raster akumulasi aliran ArcGIS --- dan saya ingin mengidentifikasi semua sel yang memiliki nilai spesifik (atau dalam kisaran nilai). Pada akhirnya, saya ingin membentuk titik yang mewakili sel-sel ini.

Saya dapat menggunakan QGIS untuk membuka hdr.adf dan mendapatkan hasil ini, alur kerjanya adalah:

• QGIS> Menu raster> Kalkulator Raster (tandai semua poin dengan nilai target)
• QGIS> Menu raster> Poligonisasi
• QGIS> Menu vektor> Submenu geometri> Centroid poligon
• Edit centroid untuk menghapus poly centroid yang tidak diinginkan (itu = 0)

Pendekatan ini "melakukan pekerjaan", tetapi tidak menarik bagi saya karena membuat 2 file yang harus saya hapus, maka saya harus menghapus catatan yang tidak diinginkan dari shapefile centroids (yaitu mereka = 0).

Sebuah pertanyaan yang ada pendekatan ini, tapi itu disesuaikan untuk ArcGIS / ArcPy, dan saya ingin tinggal di ruang FOSS.

Adakah yang memiliki resep / skrip GDAL / Python yang ada yang menginterogasi nilai sel raster, dan ketika nilai target --- atau nilai dalam rentang target --- ditemukan, catatan ditambahkan ke shapefile? Ini tidak hanya akan menghindari interaksi UI, tetapi juga akan menghasilkan hasil yang bersih dalam sekali lintasan.

Saya mengambil kesempatan itu dengan bekerja melawan salah satu presentasi Chris Garrard , tetapi pekerjaan raster tidak ada di ruang kemudi saya dan saya tidak ingin mengacaukan pertanyaan dengan kode saya yang lemah.

Jika ada yang ingin bermain dataset yang tepat, saya taruh di sini sebagai .zip .

[Edit Catatan] Meninggalkan ini untuk anak cucu. Lihat pertukaran komentar dengan om_henners. Pada dasarnya nilai x / y (baris / kolom) dibalik. Jawaban asli memiliki baris ini:

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

terbalik, seperti ini:

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

Ketika saya pertama kali menemukan masalah yang diilustrasikan dalam tangkapan layar, saya bertanya-tanya apakah saya mengimplementasikan geometri secara tidak benar, dan saya bereksperimen dengan membalik nilai koordinat x / y pada baris ini:

point.SetPoint(0, x, y)

..sebagai..

point.SetPoint(0, y, x)

Namun itu tidak berhasil. Dan saya tidak berpikir untuk mencoba membalik nilai-nilai dalam ekspresi Numpy om_henners ', percaya salah bahwa membalikkan mereka di kedua baris adalah setara. Saya pikir masalah sebenarnya berkaitan dengan nilai x_sizedan y_size, masing 30- masing dan -30, yang diterapkan ketika indeks baris dan kolom digunakan untuk menghitung titik koordinat untuk sel.

[Sunting Asli]

@om_henners, saya mencoba solusi Anda, dalam konser dengan beberapa pasangan untuk membuat shapefile titik menggunakan ogr ( invisibleroads.com , Chris Garrard ), tapi saya mengalami masalah di mana poin muncul seolah-olah tercermin melintasi garis yang lewat melalui 315/135-derajat.

Poin biru muda : dibuat oleh pendekatan QGIS saya , di atas

Poin ungu : dibuat oleh kode py GDAL / OGR , di bawah ini

[Terpecahkan]

Kode Python ini mengimplementasikan solusi lengkap seperti yang diusulkan oleh @om_henners. Saya sudah mengujinya dan berhasil. Terima kasih sobat!

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

Di GDAL Anda dapat mengimpor raster sebagai array 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)

Kemudian menggunakan numpy itu adalah masalah sederhana untuk mendapatkan indeks array yang cocok dengan ekspresi boolan:

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

Dari geotransform raster kita bisa mendapatkan informasi seperti koordinat x dan y kiri atas dan ukuran sel.

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

Sel kiri atas berhubungan dengan a[0, 0]. Ukuran Y akan selalu negatif, jadi menggunakan indeks x dan y Anda dapat menghitung koordinat setiap sel berdasarkan indeks.

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

Dari sini adalah masalah yang cukup sederhana untuk membuat shapefile menggunakan OGR. Untuk beberapa kode sampel, lihat pertanyaan ini tentang cara membuat dataset baru dengan informasi titik.

Mengapa tidak menggunakan kotak alat Sexante di QGIS? Ini seperti Model Builder untuk ArcGIS. Dengan begitu, Anda dapat rantai operasi dan memperlakukannya sebagai satu operasi. Anda dapat mengotomatiskan penghapusan file antara dan penghapusan catatan yang tidak diinginkan jika saya tidak salah.

Mungkin bermanfaat untuk mengimpor data ke postgis (dengan dukungan raster) dan menggunakan fungsi di sana. Tutorial ini mungkin memiliki elemen yang Anda butuhkan.