¿Cómo obtener las coordenadas XY y el valor de celda de cada píxel en un raster utilizando Python?

Question

¿Cómo obtener las coordenadas XY y el valor de celda de cada píxel en un raster utilizando Python?

#1 de Marcin (10 votos)
#2 de hmfly (7 votos)
#3 de Luke (7 votos)
#4 de JamesS (6 votos)
#5 de dango (3 votos)
#6 de David (3 votos)
#7 de brokev03 (2 votos)
#8 de dgketchum (2 votos)

15

Soy realmente nuevo en Python y me gustaría saber si hay un método rápido para obtener los valores de celda de un raster píxel por píxel y las coordenadas (mapa XY coordenada del centro de cada píxel) usando Python en ArcGIS 10 ?

Para describir esto más detalladamente, necesito obtener el mapa X, el mapa Y y el valor de celda del primer píxel y asignar esos tres valores a tres variables y repetir este paso para el resto de los otros píxeles (recorrer todo el ráster). ).

Creo que necesito describir mi pregunta más. El problema es que necesito obtener la ubicación X Y de un píxel del primer ráster y obtener valores de celda de varios otros rásteres correspondientes a esa ubicación X Y. Este proceso debe pasar por cada píxel del primer ráster sin crear ningún shapefile de puntos intermedios, ya que va a consumir mucho tiempo ya que tengo que manejar un ráster con casi 8 mil millones de píxeles. Además, necesito hacer esto usando Python en ArcGIS 10.

@JamesS: Muchas gracias por tu sugerencia. Sí, esto funcionaría para un ráster, pero también debo recopilar los valores de celda para otros rásteres. El problema es que, después de obtener las coordenadas X e Y del primer píxel del primer ráster, necesito obtener el valor de celda del segundo ráster correspondiente a esa ubicación X, Y del primer ráster, luego el tercer ráster y así sucesivamente. Por lo tanto, creo que cuando se recorre en bucle el primer ráster, obtener la ubicación X e Y de un píxel y obtener los valores de celda del otro ráster correspondiente a esa ubicación se debe hacer simultáneamente, pero no estoy seguro. Esto se puede hacer al convertir el primer ráster en un shapefile de puntos y realizar la función Extraer multivalores a puntos en ArcGIS 10, pero no puedo seguir ese método porque no quiero crear un shapefile ya que realmente será lento y TENGO QUE ENCONTRAR UNA SOLUCIÓN UTILIZANDO PYTHON, no con ninguna herramienta existente en ArcGIS.

@hmfly: Gracias, sí, este método (RastertoNumpyarray) funcionará si puedo obtener la coordenada de un valor de fila y columna conocido de la matriz.

@whuber: no quiero realizar ningún cálculo, todo lo que debo hacer es escribir las coordenadas XY y los valores de las celdas en un archivo de texto y eso es todo

raster python arcgis-10.0

pregunta Dash 24.02.2012 - 06:45

8 respuestas

Lea otras preguntas en las etiquetas raster python arcgis-10.0

¿Actualizando las ubicaciones de los archivos en los proyectos de QGIS? ¿Se necesita una maestría para una carrera en SIG? [cerrado]

score 10 · Answer 1

Siguiendo la idea de @ Dango, creé y probé (en pequeños rasters con la misma extensión y tamaño de celda) el siguiente código:

import arcpy, numpy

inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"

##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)

#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin

# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width

##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
    coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
    return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord

def rasCentrY(rasHeight, rasWidth):
    coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
    return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord

#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))


##Raster conversion to NumPy Array
#create NumPy array from input rasters 
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)

#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)


# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))

Según el código de @hmfly, puede tener acceso a los valores deseados:

(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
    for col in range(0,width):
        #now you have access to single array of values for one cell location

Lamentablemente hay un 'pero': el código es correcto para las matrices NumPy que pueden manejarse con la memoria del sistema. Para mi sistema (8GB), la matriz más grande fue de aproximadamente 9000,9000.

Como mi experiencia no me permite brindar más ayuda, puede considerar algunas sugerencias sobre cómo tratar con matrices grandes: enlace

El método

arcpy.RasterToNumPyArray permite especificar el subconjunto de ráster convertido a matriz NumPy ( Página de ayuda de ArcGIS10 ) lo que puede ser útil al dividir grandes conjuntos de datos en submatrices.

score 7 · Answer 2

Si solo desea obtener los valores de píxel a través de (fila, columna), puede escribir un script arcpy como este:

import arcpy
raster = arcpy.Raster("yourfilepath")
array = arcpy.RasterToNumPyArray(raster)
(height, width)=array.shape
for row in range(0,height):
    for col in range(0,width):
        print str(row)+","+str(col)+":"+str(array.item(row,col))

Pero, si desea obtener la coordenada del píxel, NumPyArray no puede ayudarlo. Puede convertir el ráster al punto mediante la herramienta RasterToPoint, y luego puede obtener la coordenada por Shape archivado.

score 7 · Answer 3

El método más sencillo para generar coordenadas y valores de celda en un archivo de texto en ArcGIS 10 es función de ejemplo , sin necesidad de código y, especialmente, no hay necesidad de recorrer cada celda. En ArcGIS < = 9.3x calculadora raster que utilizó para que sea tan simple como outfile.csv = sample(someraster) que generaría un archivo de texto de todos los valores de celda (no nulos) y coordenadas (en formato z, x, y). En ArcGIS 10, parece que el argumento "in_location_data" ahora es obligatorio, por lo que debe usar la sintaxis Sample(someraster, someraster, outcsvfile) .

Editar: También puede especificar varios rásteres: Sample([someraster, anotherraster, etc], someraster, outcsvfile) . Si esto funcionaría en 8 mil millones de células, no tengo ninguna idea ...

Editar: Nota, no he probado esto en ArcGIS 10, pero he usado la función de muestra durante años en < = 9.3 (y Estación de trabajo).

Editar: ahora he probado en ArcGIS 10 y no se enviará a un archivo de texto. La herramienta cambia la extensión del archivo a ".dbf" automáticamente. Sin embargo ... el siguiente código de Python funciona como las declaraciones de álgebra de mapas aún son compatibles con ArcGIS 10:

import arcgisscripting
gp=arcgisscripting.create()
gp.multioutputmapalgebra(r'%s=sample(%s)' % (outputcsv,inputraster))

score 6 · Answer 4

Una forma de hacer esto sería utilizar Raster_To_Point herramienta seguida de la herramienta Add_XY_Coordinates . Terminará con un shapefile donde cada fila en la tabla de atributos representa un píxel de su ráster con columnas para X_Coord , Y_Coord y Cell_Value . Luego puede recorrer esta tabla con un cursor (o exportarlo a algo como Excel, si lo prefiere).

Si solo tienes un ráster para procesar, probablemente no valga la pena hacer un script, solo usa las herramientas de ArcToolbox. Si necesita hacer esto para muchos rásteres, podría intentar algo como esto:

[ Nota: No tengo ArcGIS 10 y no estoy familiarizado con ArcPy, por lo que este es solo un esquema muy aproximado. No se ha probado y es casi seguro que se necesitarán ajustes para que funcione.]

import arcpy, os
from arcpy import env

# User input
ras_fold = r'path/to/my/data'           # The folder containing the rasters
out_fold = r'path/to/output/shapefiles' # The folder in which to create the shapefiles

# Set the workspace
env.workspace = ras_fold

# Get a list of raster datasets in the raster folder
raster_list = arcpy.ListRasters("*", "All")

# Loop over the rasters
for raster in raster_list:
    # Get the name of the raster dataset without the file extension
    dataset_name = os.path.splitext(raster)[0]

    # Build a path for the output shapefile
    shp_path = os.path.join(out_fold, '%s.shp' % dataset_name)

    # Convert the raster to a point shapefile
    arcpy.RasterToPoint_conversion(raster, shp_path, "VALUE")

    # Add columns to the shapefile containing the X and Y co-ordinates
    arcpy.AddXY_management(shp_path)

Luego puede recorrer las tablas de atributos de shapefile usando un Cursor de búsqueda o (posiblemente más sencillo) utilizando dbfpy . Esto le permitirá leer los datos de su ráster (ahora almacenados en una tabla .dbf de shapefile) en las variables de python.

from dbfpy import dbf

# Path to shapefile .dbf
dbf_path = r'path\to\my\dbf_file.dbf'

# Open the dbf file
db = dbf.Dbf(dbf_path)

# Loop over the records
for rec in db:
    cell_no = rec['POINTID'] # Numbered from top left, running left to right along each row
    cell_x = rec['POINT_X']
    cell_y = rec['POINT_Y']
    cell_val = rec['GRID_CODE']

    # Print values
    print cell_no, cell_x, cell_y, cell_val

score 3 · Answer 5

Tal vez podría crear un archivo de mundo para el ráster, convertir el ráster en una matriz numpy. Luego, si recorres la matriz, obtendrás los valores de las celdas y si de forma increamentaria, actualice x, y del archivo mundial, también tendrá las coordenadas para cada valor de celda. Espero que sea de utilidad.

score 3 · Answer 6

El código de Marcin funcionó bien, excepto que un problema en las funciones rasCentrX y rasCentrY hacía que las coordenadas de salida aparecieran en una resolución diferente (como observó Grazia). Mi solución fue cambiar

coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)

a

coordX = rasXmin + ((0.5 + rasWidth) * rasMeanCellWidth)

y

  coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)

a

  coordY = rasYmin + ((0.5 + rasHeight) * rasMeanCellHeight)

Utilicé el código para convertir una cuadrícula ESRI en un archivo CSV. Esto se logró eliminando la referencia a inRaster2, luego usando un csv.writer para generar las coordenadas y los valores:

out = csv.writer(open(outputCSV,"wb"), delimiter=',', quoting=csv.QUOTE_NONNUMERIC)
out.writerow(['X','Y','Value'])
(height, width, dim )=inRasterFullArray.shape
for row in range(0,height):
    for col in range(0,width):
        out.writerow(inRasterFullArray[row,col])

Tampoco encontré la transposición que se necesitaba en

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

así lo convirtió a

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray))

score 2 · Answer 7

Feo pero altamente efectivo:

Cree una nueva entidad de puntos con 4 puntos fuera de las esquinas del ráster en cuestión. Asegúrese de utilizar el mismo sistema de coordenadas que el ráster en cuestión.
Agregue los campos dobles 'xcor' y 'ycor'
Calcular geometría para obtener coordenadas para estos campos
Spatial Analyst- > Interpolation- > Trend - > Regresión lineal
Configuración del entorno: ajustar el ráster y el tamaño de celda al mismo que el ráster en cuestión
Ejecutar por separado para 'xcor' y 'ycor'
Los clasificadores de salida vienen con coordenadas como valores de celda, se usan como entrada para scripts.

score 2 · Answer 8

Una solución simple que utiliza paquetes de código abierto de python:

import fiona
import rasterio
from pprint import pprint


def raster_point_coords(raster, points):

    # initialize dict to hold data
    pt_data = {}

    with fiona.open(points, 'r') as src:
        for feature in src:
            # create dict entry for each feature
            pt_data[feature['id']] = feature

    with rasterio.open(raster, 'r') as src:
        # read raster into numpy array
        arr = src.read()
        # rasterio always reads into 3d array, this is 2d, so reshape
        arr = arr.reshape(arr.shape[1], arr.shape[2])
        # get affine, i.e. data needed to work between 'image' and 'raster' coords
        a = src.affine

    for key, val in pt_data.items():
        # get coordinates
        x, y = val['geometry']['coordinates'][0], val['geometry']['coordinates'][1]
        # use affine to convert to row, column
        col, row = ~a * (x, y)
        # remember numpy array is indexed array[row, column] ie. y, x
        val['raster_value'] = arr[int(row), int(col)]

    pprint(pt_data) 

if __name__ == '__main__':
    # my Landsat raster
    ras = '/data01/images/sandbox/LT05_040028_B1.tif'
    # my shapefile with two points which overlap raster area
    pts = '/data01/images/sandbox/points.shp'
    # call function
    raster_point_coords(ras, pts)

Fiona es útil ya que puedes abrir un shapefile, recorrer las funciones y (como yo lo he hecho) agregarlas a un objeto dict . De hecho, el propio Fiona feature también es como un dict , por lo que es fácil acceder a las propiedades. Si mis puntos tuvieran algún atributo, aparecerían en este dictado junto con las coordenadas, id, etc.

Rasterio es útil porque es fácil de leer en la trama como una matriz numpy, un tipo de datos ligero y rápido. También tenemos acceso a un dict de las propiedades del ráster, incluido el affine , que es todos los datos que necesitamos para convertir las coordenadas del ráster x, y en la fila de la matriz, coordenadas col. Consulte la excelente explicación de @ perrygeo aquí .

Terminamos con un pt_data del tipo dict que tiene datos para cada punto y el raster_value extraído. Podríamos reescribir fácilmente el shapefile con los datos extraídos también si quisiéramos.