GDAL y Python: ¿Cómo obtener las coordenadas de todas las celdas que tienen un valor específico?

Question

Tengo un Rejilla binaria Arc/Info ---específicamente, un raster de acumulación de flujo de ArcGIS--- y me gustaría identificar todas las celdas que tienen un valor específico (o en un rango de valores). En última instancia, me gustaría un shapefile de puntos que representan estas células.

Puedo utilizar QGIS para abrir el hdr.adf y obtener este resultado, el flujo de trabajo es:

• QGIS > Menú Raster > Calculadora Raster (marcar todos los puntos con el valor objetivo)
• QGIS > Menú Raster > Poligonizar
• QGIS > Menú Vectorial > Submenú Geometría > Centros de polígonos
• Editar los centroides para eliminar los poli centros no deseados (aquellos = 0)

Este enfoque "hace el trabajo", pero no me atrae porque crea 2 archivos que tengo que borrar, y luego tengo que eliminar los registros no deseados del shapefile de centroides (es decir, aquellos = 0).

Un pregunta existente aborda este tema, pero está adaptado a ArcGIS/ArcPy, y me gustaría permanecer en el espacio FOSS.

¿Alguien tiene una receta/script de GDAL/Python que interrogue los valores de las celdas de un raster, y cuando se encuentre un valor objetivo -o un valor en un rango objetivo- se añada un registro a un shapefile? Esto no sólo evitaría la interacción con la interfaz de usuario, sino que crearía un resultado limpio en una sola pasada.

Lo intenté por medio de trabajando contra una de las presentaciones de Chris Garrard pero el trabajo de rasterización no es mi especialidad y no quiero saturar la pregunta con mi débil código.

Si alguien quiere el conjunto de datos exacto para jugar, Lo pongo aquí como un .zip .

---

[Editar Notas] Dejando esto para la posteridad. Ver los intercambios de comentarios con om_henners. Básicamente, los valores x/y (fila/columna) estaban invertidos. La respuesta original tenía esta línea:

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

invertido, así:

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

Cuando encontré por primera vez el problema ilustrado en la captura de pantalla, me pregunté si había implementado la geometría de forma incorrecta, y experimenté invirtiendo los valores de las coordenadas x/y en esta línea:

point.SetPoint(0, x, y)

..como..

point.SetPoint(0, y, x)

Sin embargo, eso no funcionó. Y no se me ocurrió intentar voltear los valores en la expresión Numpy de om_henners, creyendo erróneamente que voltearlos en cualquier línea era equivalente. Creo que el verdadero problema está relacionado con el x_size y y_size valores, respectivamente 30 y -30 que se aplican cuando los índices de fila y columna se utilizan para calcular las coordenadas de los puntos de las celdas.

[Edición original]

@om_henners, estoy probando tu solución, en conjunto con un par de recetas para hacer shapefiles de puntos usando ogr ( invisibleroads.com , Chris Garrard ), pero tengo un problema en el que los puntos aparecen como si se reflejaran en una línea que pasa por 315/135 grados.

Puntos azul claro : creado por mi Enfoque de QGIS , arriba

Puntos morados : creado por el Código GDAL/OGR py , a continuación

---

[Resuelto]

Este código de Python implementa la solución completa propuesta por @om_henners. Lo he probado y funciona. ¡Gracias tío!

---

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

Answer 1

En GDAL puedes importar el raster como un 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)

Entonces, utilizando numpy es muy sencillo obtener los índices de un array que coincida con una expresión boolan:

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

De la geotransformación del raster podemos obtener información como las coordenadas x e y superiores izquierdas y los tamaños de las celdas.

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

La celda superior izquierda corresponde a a[0, 0] . El tamaño Y siempre será negativo, así que usando los índices x e y puedes calcular las coordenadas de cada celda basándote en los í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 de aquí es bastante sencillo crear un shapefile utilizando OGR. Para ver algunos ejemplos de código ver esta pregunta para saber cómo generar un nuevo conjunto de datos con información puntual.

Answer 2

¿Por qué no utilizar el Caja de herramientas Sexante en QGIS? Es algo así como el Model Builder de ArcGIS. De esta manera, usted puede encadenar las operaciones y tratarlas como una sola operación. Puedes automatizar la eliminación de archivos intermedios y la eliminación de registros no deseados si no me equivoco.

Answer 3

Puede ser útil importar los datos a Postgis (con soporte raster) y utilizar las funciones allí. Este tutorial puede tener elementos que necesitas.