¿Cómo iterar a través de cada célula en un continuo de trama?

Question

Consulte este enlace para obtener más detalles.

El Problema:

Quiero un bucle a través de un continuo de trama (uno que no tiene tabla de atributos), celda por celda, y obtener el valor de la celda. Quiero tomar esos valores y ejecutar las oraciones condicionales en ellos, emulando el mapa de álgebra pasos que se detallan a continuación sin que en realidad el uso de la calculadora ráster.

Por petición de comentarios de abajo, he añadido detalles proporcionar antecedentes del problema y justificación de la necesidad de implementar un método como tal en la sección denominada "El análisis necesario:".

El análisis que se propone a continuación, siendo relevante para mi problema por proporcionar el fondo,no necesita ser implementado en una respuesta. El alcance de la pregunta, se refiere a la iteración de un continuo de trama para obtener o establecer los valores de celda.

El análisis necesario:

Si ALGUNA de las siguientes condiciones se cumplen, dar celda de salida a un valor de 1. Sólo dar celda de salida a un valor de 0 si ninguna de las condiciones se cumplen.

Condición 1: Si el valor de la celda es mayor que la parte superior e inferior de las células, dar valor de 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Donde archivo del kernel se parece a esto:

3 3 
0 1 0
0 0 0
0 1 0

Condición 2: Si el valor de la celda es mayor que la de la izquierda y la derecha de las células, dar valor de 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Donde archivo del kernel se parece a esto:

3 3 
0 0 0
1 0 1
0 0 0  

Condición 3: Si el valor de la celda es mayor que topleft y bottomright células, dar valor de 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Donde archivo del kernel se parece a esto:

3 3 
1 0 0
0 0 0
0 0 1 

Condición 4: Si el valor de la celda es mayor que bottomleft y topright células, dar valor de 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Donde archivo del kernel se parece a esto:

3 3 
0 0 1
0 0 0
1 0 0 

Condición 5: Si cualquier una de las células adyacentes que tiene un valor IGUAL a el centro de la célula, dar el ráster de salida un valor de 1 (uso de focales variedad con dos más cercano barrio de los cálculos)

¿Por qué no utilizar el mapa de álgebra?

Se ha señalado por debajo de la que mi problema puede ser resuelto mediante la asignación de álgebra, pero como se ha visto anteriormente este es un gran total de seis trama cálculos, más de uno de combinar todos los rásteres creado juntos. A mí me parece que es mucho más eficiente para ir celda por celda y hacer todas las comparaciones a la vez en cada celda en lugar de un bucle a través de cada uno de forma individual siete veces y utilizando un poco de memoria para crear siete rásteres.

¿Cómo debe el problema de ser atacado?

El enlace de arriba aconseja el uso de IPixelBlock interfaz, sin embargo no está claro a partir de ESRI documentación que usted tenga acceso a un único valor de la celda en sí a través de IPixelBlock, o si tiene acceso a varios valores de celda del tamaño de la IPixelBlock establecido. Una buena respuesta debe sugerir un método para acceder a los valores de celda de un continuo de trama y de proporcionar una explicación de la metodología detrás de el código, si no lo aparentemente obvio.

En resumen:

¿Cuál es el mejor método de bucle a través de cada célula en un CONTINUO de trama (el cual tiene ningún atributo de la tabla) para acceder a sus valores de celda?

Una buena respuesta, no necesita implementar el análisis de los pasos descritos anteriormente, se necesita sólo para proporcionar una metodología para acceder a los valores de las celdas de una trama.

Answer 1

Veo que esto ya ha sido resuelto por el Cartel Original (OP), pero voy a publicar una sencilla solución en python sólo en caso de que alguien en el futuro está interesado en las diferentes formas de resolver este problema. Yo soy parcial a software de código abierto, así que aquí es una solución de uso de GDAL en python:

import gdal

#Set GeoTiff driver
driver = gdal.GetDriverByName("GTiff")
driver.Register()

#Open raster and read number of rows, columns, bands
dataset = gdal.Open(filepath)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
allBands = dataset.RasterCount
band = dataset.GetRasterBand(1)

#Get array of raster cell values.  The two zeros tell the 
#iterator which cell to start on and the 'cols' and 'rows' 
#tell the iterator to iterate through all columns and all rows.
def get_raster_cells(band,cols,rows):
    return band.ReadAsArray(0,0,cols,rows)

Implementar la función como esta:

#Bind array to a variable
rasterData = get_raster_cells(band,cols,rows)

#The array will look something like this if you print it
print rasterData
> [[ 1, 2, 3 ],
   [ 4, 5, 6 ],
   [ 7, 8, 9 ]]

A continuación, recorrer sus datos con un bucle anidado:

for row in rasterData:
    for val in row:
        print val
> 1
  2
  3
  4...

O tal vez usted quiere acoplar el 2-D de la matriz con una lista de comprensión:

flat = [val for row in rasterData for val in row]

De todos modos, mientras la iteración a través de los datos en una celda por celda es posible arrojar algo de condicionales en su bucle para cambiar/editar los valores. Ver este guión que escribió para diferentes formas de acceso a los datos: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py.

Answer 2

Actualización! El numpy solución:

import arcpy
import numpy as np

in_ras = path + "/rastername"

raster_Array = arcpy.RasterToNumPyArray(in_ras)
row_num = raster_Array.shape[0]
col_num = raster_Array.shape[1]
cell_count = row_num * row_num

row = 0
col = 0
temp_it = 0

while temp_it < cell_count:
    # Insert conditional statements
    if raster_Array[row, col] > 0:
        # Do something
        val = raster_Array[row, col]
        print val
    row+=1
    if col > col_num - 1:
        row = 0
        col+=1

Así, conseguir que el acabado de la matriz de vuelta a ráster utilizando arcpy es problemático. arcpy.NumPyArrayToRaster es squirrelly y tiende a redefinir las extensiones incluso si usted alimenta a su LL coordenadas.

Prefiero guardar como texto.

np.savetxt(path + "output.txt", output, fmt='%.10f', delimiter = " ")

Estoy ejecutando Python como de 64 bits para la velocidad - a partir de ahora esto significa que no puede alimentar a numpy.savetxt un encabezado. Así que tengo que abrir la salida y agregar el ASCII de encabezado que el Arco quiere que antes de la conversión de ASCII a Trama

File_header = "NCOLS xxx" + '\n'+ "NROWS xxx" + '\n' + "XLLCORNER xxx"+'\n'+"YLLCORNER xxx"+'\n'+"CELLSIZE xxx"+'\n'+"NODATA_VALUE xxx"+'\n'

El numpy versión corre mi cambio de trama, multiplicaciones, y además mucho más rápido (1000 iteraciones en 2 minutos) de la arcpy versión (1000 iteraciones en 15 min)

VERSIÓN ANTIGUA me puede eliminar esta tarde Acabo de escribir un guión similar. He tratado de convertir a los puntos y a través de la búsqueda de cursor. Tengo sólo 5000 iteraciones en 12 horas. Busqué, entonces, de otra manera.

Mi manera de hacer esto es recorrer el centro de la celda de coordenadas de cada celda. Me inicio en la esquina superior izquierda y mover de derecha a izquierda. Al final de la fila I mover hacia abajo una fila y empezar de nuevo desde la izquierda. Tengo un 240 m ráster con 2603 columnas y 2438 filas para un total de 6111844 total de células. Yo uso una variable de iteración y de un bucle while. Ver más abajo

Un par de notas: 1 - es necesario conocer las coordenadas de la medida

2 - ejecutar con el punto de las coordenadas de centro de la celda - se mueven en 1/2 el tamaño de la celda a partir de la medida de los valores de

3 - Mi script utiliza el valor de la celda para tirar de un valor específico de la trama, a continuación, mayús esta trama se centre en la célula original. Esto añade un cero trama para ampliar el alcance antes de añadir en un final de trama. Este es sólo un ejemplo. Usted puede poner sus instrucciones condicionales en aquí (segundo si la instrucción dentro del bucle while).

4 - Este script asume todos los de la trama valores pueden ser emitidos como enteros. Esto significa que usted necesita para deshacerse de los datos en primer lugar. Con IsNull.

6 - todavía no estoy feliz con esto y estoy trabajando para tomar esta fuera de arcpy completamente. Prefiero el reparto de la arrays de numpy y hacer las matemáticas más allá luego traerlo de vuelta a Arco.

ULx = 959415 ## coordinates for the Upper Left of the entire raster 
ULy = 2044545
x = ULx ## I redefine these if I want to run over a smaller area
y = ULy
temp_it = 0

while temp_it < 6111844: # Total cell count in the data extent
        if x <= 1583895 and y >= 1459474: # Coordinates for the lower right corner of the raster
           # Get the Cell Value
           val_result = arcpy.GetCellValue_management(inraster, str(x)+" " +str(y), "1")
           val = int(val_result.getOutput(0))
        if val > 0: ## Here you could insert your conditional statements
            val_pdf = Raster(path + "pdf_"str(val))
            shift_x  =  ULx - x # This will be a negative value
            shift_y = ULy - y # This will be a positive value
            arcpy.Shift_management(val_pdf, path+ "val_pdf_shift", str(-shift_x), str(-shift_y))
            val_pdf_shift = Raster(path + "val_pdf_shift")
            val_pdf_sh_exp = CellStatistics([zeros, val_pdf_shift], "SUM", "DATA")
            distr_days = Plus(val_pdf_sh_exp, distr_days)
        if temp_it % 20000 == 0: # Just a print statement to tell me how it's going
                print "Iteration number " + str(temp_it) +" completed at " + str(time_it)
        x += 240 # shift x over one column
        if x > 1538295: # if your at the right hand side of a row
            y = y-240 # Shift y down a row
            x = 959415 # Shift x back to the first left hand column
        temp_it+=1

distr_days.save(path + "Final_distr_days")

Answer 3

Trate de usar IGridTable, ICursor, IRow. Este fragmento de código es para la actualización de los valores de celda ráster, sin embargo, se muestran los conceptos básicos de la iteración:

¿Cómo puedo agregar un nuevo campo en una tabla de atributos de ráster y bucle a través de él?

Public Sub CalculateArea(raster As IRaster, areaField As String)
    Dim bandCol As IRasterBandCollection
    Dim band As IRasterBand

    Set bandCol = raster
    Set band = bandCol.Item(0)

    Dim hasTable As Boolean
    band.hasTable hasTable
    If (hasTable = False) Then
        Exit Sub
    End If    

    If (AddVatField(raster, areaField, esriFieldTypeDouble, 38) = True) Then
        ' calculate cell size
        Dim rstProps As IRasterProps
        Set rstProps = raster

        Dim pnt As IPnt
        Set pnt = rstProps.MeanCellSize

        Dim cellSize As Double
        cellSize = (pnt.X + pnt.Y) / 2#

        ' get fields index
        Dim attTable As ITable
        Set attTable = band.AttributeTable

        Dim idxArea As Long, idxCount As Long
        idxArea = attTable.FindField(areaField)
        idxCount = attTable.FindField("COUNT")

        ' using update cursor
        Dim gridTableOp As IGridTableOp
        Set gridTableOp = New gridTableOp

        Dim cellCount As Long, cellArea As Double

        Dim updateCursor As ICursor, updateRow As IRow
        Set updateCursor = gridTableOp.Update(band.RasterDataset, Nothing, False)
        Set updateRow = updateCursor.NextRow()
        Do Until updateRow Is Nothing
            cellCount = CLng(updateRow.Value(idxCount))
            cellArea = cellCount * (cellSize * cellSize)

            updateRow.Value(idxArea) = cellArea
            updateCursor.updateRow updateRow

            Set updateRow = updateCursor.NextRow()
        Loop

    End If
End Sub

Una vez que están recorriendo la tabla se puede obtener el campo específico de la fila de valor mediante el uso de row.get_Value(yourfieldIndex). Si buscas en Google

arcobjects fila.get_Value

usted debe ser capaz de conseguir un montón de ejemplos que muestran esto.

Espero que ayude.

Answer 4

¿Y esto como una idea radical, que se requieren para programar en python o de ArcObjects.

1. Convertir su grilla a un punto featureclasss.
2. Crear XY campos y rellenar.
3. La carga de los puntos en un diccionario donde la clave es una cadena de X,Y y elemento es el valor de la celda..
4. El paso a través de su diccionario, y para cada punto de trabajo el 8 rodean de células Xy.
5. Recuperar estas de su diccionario y test con sus reglas, tan pronto como se encuentre un valor que es cierto que usted puede saltarse el resto de las pruebas.
6. Escribir los resultados en otro diccionario y, a continuación, volver a convertir a una cuadrícula para crear un punto de FeatureClass y, a continuación, convertir puntos en una cuadrícula.

Answer 5

AFAIK datos ráster se puede leer en tres formas:

• por celda (ineficaces);
• por imagen (muy eficiente);
• por bloques (de la forma más eficiente).

Sin reinventar la rueda, sugiero leer estas iluminadoras las diapositivas de Chris Garrard.

Así que el método más eficiente es leer los datos por bloque, sin embargo, esto podría causar una pérdida de datos en la correspondencia de los píxeles ubicados sobre el bloque de límites al aplicar el filtro. Por lo que una alternativa segura forma en que debe consistir en la lectura de la imagen completa de una sola vez y mediante el numpy enfoque.

En el cómputo de lado, en cambio, que debo usar gdalfilter.py e implícitamente la VRT KernelFilteredSource enfoque en el fin de aplicar los filtros necesarios y, sobre todo, evitar pesados cálculos.