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