GDAL e Python: como obter coordenadas para todas as células com um valor específico?

Eu tenho uma grade binária Arc / Info - especificamente, uma varredura de acumulação de fluxo ArcGIS - e gostaria de identificar todas as células com um valor específico (ou em um intervalo de valores). Por fim, eu gostaria de um shapefile de pontos representando essas células.

Posso usar o QGIS para abrir o hdr.adf e obter esse resultado, o fluxo de trabalho é:

• QGIS> menu Raster> Calculadora Raster (marque todos os pontos com o valor alvo)
• QGIS> menu Raster> poligonizar
• QGIS> Menu vetorial> Submenu Geometria> Centróide de polígono
• Edite os centróides para excluir os poli-centróides indesejados (aqueles = 0)

Essa abordagem "faz o trabalho", mas não me agrada porque cria 2 arquivos que tenho que excluir e, em seguida, tenho que remover registros indesejados do shapefile dos centróides (ou seja, aqueles = 0).

Uma pergunta existente aborda esse assunto, mas é adaptada ao ArcGIS / ArcPy, e eu gostaria de permanecer no espaço do FOSS.

Alguém tem uma receita / script GDAL / Python existente que interroga os valores das células de uma varredura e, quando um valor alvo --- ou um valor no intervalo alvo --- é encontrado, um registro é adicionado a um shapefile? Isso não apenas evitaria a interação da interface do usuário, mas também criaria um resultado limpo em uma única passagem.

Tomei uma chance trabalhando contra uma das apresentações de Chris Garrard , mas o trabalho raster não está na minha casa do leme e não quero bagunçar a questão com meu código fraco.

Se alguém quiser jogar com o conjunto de dados exato, coloquei-o aqui como um .zip .

[Editar notas] Deixando isso para trás para a posteridade. Veja trocas de comentários com om_henners. Basicamente, os valores x / y (linha / coluna) foram invertidos. A resposta original tinha esta linha:

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

invertido, assim:

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

Quando encontrei o problema ilustrado na captura de tela, imaginei se havia implementado a geometria incorretamente e experimentei inverter os valores das coordenadas x / y nesta linha:

point.SetPoint(0, x, y)

..Como..

point.SetPoint(0, y, x)

No entanto, isso não funcionou. E não pensei em tentar inverter os valores na expressão Numpy do om_henners, acreditando erroneamente que lançá-los em qualquer linha era equivalente. Penso que o verdadeiro problema está relacionado com os valores x_sizee y_size, respectivamente, 30e -30que são aplicados quando os índices de linha e coluna são usados ​​para calcular coordenadas de ponto para as células.

[Edição original]

@om_henners, estou tentando sua solução, em conjunto com algumas receitas para criar arquivos shapefiles de pontos usando ogr ( invisibleroads.com , Chris Garrard ), mas estou tendo um problema em que os pontos aparecem como se refletidos em uma linha que passava a 315/135 graus.

Pontos azuis claros : criados pela minha abordagem QGIS , acima

Pontos roxos : criados pelo código py GDAL / OGR , abaixo

[Resolvido]

Este código Python implementa a solução completa, conforme proposto por @om_henners. Eu testei e funciona. Valeu cara!

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

No GDAL, você pode importar a varredura como uma matriz 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)

Então, usando numpy, é simples obter os índices de uma matriz que correspondem a uma expressão boolan:

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

A partir da geotransformação raster, podemos obter informações como as coordenadas xey superiores e os tamanhos das células.

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

A célula superior esquerda corresponde a a[0, 0]. O tamanho Y sempre será negativo; portanto, usando os índices x e y, é possível calcular as coordenadas de cada célula com base nos índices.

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

A partir daqui, é bastante simples criar um shapefile usando OGR. Para algum código de exemplo, consulte esta pergunta para saber como gerar um novo conjunto de dados com informações de ponto.

Por que não usar a caixa de ferramentas Sexante no QGIS? É como o Model Builder for ArcGIS. Dessa forma, você pode encadear operações e tratá-las como uma operação. Você pode automatizar a exclusão de arquivos intermediários e a remoção de registros indesejados, se não me engano.

Pode ser útil importar os dados para o postgis (com suporte a varredura) e usar as funções lá. Este tutorial pode ter elementos que você precisa.