14 votos

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

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

enter image description here


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

2 votos

Un consejo rápido para su código: puede utilizar la proyección rasterizada como su proyección shapefile con srs.ImportFromWkt(r.GetProjection()) (en lugar de tener que crear una proyección a partir de una cadena proj conocida).

0 votos

Tu código hace todo lo que busco, excepto que no incluye el valor de la celda raster, ya que has escrito tu filtro numpy para que sólo incluya valores = 1000. ¿Podría indicarme la dirección correcta para incluir los valores de las celdas rasterizadas en la salida? Gracias.

23voto

ESV Puntos 4591

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.

0 votos

Hola amigo, estoy teniendo un pequeño problema para implementar esto. He actualizado la pregunta para incluir el código que estoy utilizando y una captura de pantalla de lo que estoy recibiendo. Básicamente, el código .py está creando una imagen especular (colocación de puntos) de lo que el enfoque de QGIS está generando. Los puntos en mi implementación caen fuera de los límites de la trama, por lo que el problema tiene que ser con mi código. := Espero que pueda arrojar algo de luz. Gracias.

0 votos

Lo siento, ha sido culpa mía. Cuando se importa un raster en GDAL las filas son la dirección y y las columnas son la dirección x. He actualizado el código de arriba, pero el truco es obtener los índices con (y_index, x_index) = np.nonzero(a > threshold)

1 votos

Además, por si acaso, observe la adición de la mitad del tamaño de la celda a las coordenadas en ambas direcciones para centrar el punto en la celda.

4voto

shsteimer Puntos 8749

¿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.

0voto

JaakL Puntos 786

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

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X