¿Cómo itero a través de cada celda en un ráster continuo?

13

Consulte este enlace para obtener más detalles.

El problema:

Quiero recorrer un continuo raster (uno que no tiene tabla de atributos ), celda por celda, y obtener el valor de la celda. Quiero tomar esos valores y ejecutar condicionales en ellos, emulando los pasos de álgebra del mapa que se detallan a continuación sin usar la calculadora raster.

A solicitud de comentarios a continuación, he agregado detalles que brindan antecedentes del problema y justifican la necesidad de implementar un método como tal en la siguiente sección llamada "El análisis necesario:".

El análisis propuesto a continuación, si bien es relevante para mi problema al proporcionar antecedentes, no necesita implementarse en una respuesta. El alcance de la pregunta solo se refiere a la iteración a través de un ráster continuo para obtener / Establecer los valores de celda.

El análisis necesario:

Si se cumple CUALQUIERA de las siguientes condiciones, asigne a la celda de salida un valor de 1. Solo dé a la celda de salida un valor de 0 si no cumple ninguna de las condiciones.

Condición 1: si el valor de la celda es mayor que las celdas superior e inferior, indique el valor 1:

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

Donde el archivo del kernel se ve así:

3 3 
0 1 0
0 0 0
0 1 0

Condición 2: si el valor de la celda es mayor que las celdas izquierda y derecha, indique el valor 1:

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

Donde el archivo del kernel se ve así:

3 3 
0 0 0
1 0 1
0 0 0  

Condición 3: si el valor de la celda es mayor que el de las celdas de fondo y de fondo, indique el valor de 1:

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

Donde el archivo del kernel se ve así:

3 3 
1 0 0
0 0 0
0 0 1 

Condición 4: Si el valor de la celda es mayor que el de las celdas de fondo y en el derecho, indique el valor de 1:

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

Donde el archivo del kernel se ve así:

3 3 
0 0 1
0 0 0
1 0 0 

Condición 5: si cualquiera una de las celdas adyacentes tiene un valor IGUAL a la celda central, asigne a la trama de salida un valor de 1 ( utilizando una variedad focal con dos cálculos de vecindarios más cercanos )

¿Por qué no usar el álgebra de mapas?

Se ha señalado a continuación que mi problema podría resolverse utilizando el álgebra de mapas, pero como se ve más arriba, se trata de un total de seis cálculos ráster, más uno para combinar todos los rásteres creados juntos. Me parece que es mucho más eficiente ir celda por celda y hacer todas las comparaciones a la vez en cada celda en lugar de recorrer cada una de ellas individualmente siete veces y utilizar bastante memoria para crear siete rásteres.

¿Cómo se debe atacar el problema?

El enlace anterior aconseja utilizar la interfaz IPixelBlock, sin embargo, no está claro en la documentación de ESRI si realmente está accediendo a un solo valor de celda a través de IPixelBlock, o si está accediendo a varios valores de celda del tamaño del IPixelBlock que estableció. Una buena respuesta debería sugerir un método para acceder a los valores de celda de un ráster continuo y proporcionar una explicación de la metodología detrás del código, si no es aparentemente obvia.

En resumen:

¿Cuál es el mejor método para recorrer cada celda en un ráster CONTINUO (que no tiene tabla de atributos ) para acceder a sus valores de celda?

Una buena respuesta no necesita implementar los pasos de análisis descritos anteriormente, solo necesita proporcionar una metodología para acceder a los valores de celda de un ráster.

    
pregunta Conor 13.04.2017 - 14:33

6 respuestas

2

Una solución:

Resolví esto antes hoy. El código es una adaptación de este método . El concepto detrás de esto no fue terriblemente difícil una vez que descubrí lo que realmente hacen los objetos utilizados para interactuar con la trama. El método siguiente toma dos conjuntos de datos de entrada (inRasterDS y outRasterDS). Ambos son el mismo conjunto de datos, acabo de hacer una copia de inRasterDS y la pasé al método como outRasterDS. De esta manera, ambos tienen la misma extensión, referencia espacial, etc. El método lee los valores de inRasterDS, celda por celda, y realiza comparaciones de vecinos más cercanos en ellos. Utiliza los resultados de esas comparaciones como valores almacenados en outRasterDS.

El proceso:

Utilicé IRasterCursor - > IPixelBlock - > SafeArray para obtener los valores de píxel e IRasterEdit para escribir nuevos valores en el ráster. Cuando crea IPixelBlock, le está diciendo a la máquina el tamaño y la ubicación del área en la que desea leer / escribir. Si solo desea editar la mitad inferior de un ráster, establezca eso como sus parámetros de IPixelBlock. Si desea recorrer todo el ráster, debe establecer IPixelBlock igual al tamaño de todo el ráster. Hago esto en el método a continuación pasando el tamaño a IRasterCursor (pSize) y luego obtengo el PixelBlock desde el cursor de la trama.

La otra clave es que tiene que usar SafeArray para interactuar con los valores de este método. Obtiene IPixelBlock de IRasterCursor, luego SafeArray de IPixelBlock. Luego lees y escribes en SafeArray. Cuando termine de leer / escribir en SafeArray, vuelva a escribir su SafeArray completo en IPixelBlock, luego escriba su IPixelBlock en IRasterCursor, luego finalmente use IRasterCursor para configurar la ubicación para comenzar a escribir e IRasterEdit para hacer la escritura. Este último paso es donde realmente editas los valores del conjunto de datos.

    public static void CreateBoundaryRaster(IRasterDataset2 inRasterDS, IRasterDataset2 outRasterDS)
    {
        try
        {
            //Create a raster. 
            IRaster2 inRaster = inRasterDS.CreateFullRaster() as IRaster2; //Create dataset from input raster
            IRaster2 outRaster = outRasterDS.CreateFullRaster() as IRaster2; //Create dataset from output raster
            IRasterProps pInRasterProps = (IRasterProps)inRaster;
            //Create a raster cursor with a pixel block size matching the extent of the input raster
            IPnt pSize = new DblPnt();
            pSize.SetCoords(pInRasterProps.Width, pInRasterProps.Height); //Give the size of the raster as a IPnt to pass to IRasterCursor
            IRasterCursor inrasterCursor = inRaster.CreateCursorEx(pSize); //Create IRasterCursor to parse input raster 
            IRasterCursor outRasterCursor = outRaster.CreateCursorEx(pSize); //Create IRasterCursor to parse output raster
            //Declare IRasterEdit, used to write the new values to raster
            IRasterEdit rasterEdit = outRaster as IRasterEdit;
            IRasterBandCollection inbands = inRasterDS as IRasterBandCollection;//set input raster as IRasterBandCollection
            IRasterBandCollection outbands = outRasterDS as IRasterBandCollection;//set output raster as IRasterBandCollection
            IPixelBlock3 inpixelblock3 = null; //declare input raster IPixelBlock
            IPixelBlock3 outpixelblock3 = null; //declare output raster IPixelBlock
            long blockwidth = 0; //store # of columns of raster
            long blockheight = 0; //store # of rows of raster

            //create system array for input/output raster. System array is used to interface with values directly. It is a grid that overlays your IPixelBlock which in turn overlays your raster.
            System.Array inpixels; 
            System.Array outpixels; 
            IPnt tlc = null; //set the top left corner

            // define the 3x3 neighborhood objects
            object center;
            object topleft;
            object topmiddle;
            object topright;
            object middleleft;
            object middleright;
            object bottomleft;
            object bottommiddle;
            object bottomright;

            long bandCount = outbands.Count; //use for multiple bands (only one in this case)

            do
            {

                inpixelblock3 = inrasterCursor.PixelBlock as IPixelBlock3; //get the pixel block from raster cursor
                outpixelblock3 = outRasterCursor.PixelBlock as IPixelBlock3;
                blockwidth = inpixelblock3.Width; //set the # of columns in raster
                blockheight = inpixelblock3.Height; //set the # of rows in raster
                outpixelblock3.Mask(255); //set any NoData values

                for (int k = 0; k < bandCount; k++) //for every band in raster (will always be 1 in this case)
                {
                    //Get the pixel array.
                    inpixels = (System.Array)inpixelblock3.get_PixelData(k); //store the raster values in a System Array to read
                    outpixels = (System.Array)outpixelblock3.get_PixelData(k); //store the raster values in a System Array to write
                    for (long i = 1; i < blockwidth - 1; i++) //for every column (except outside columns)
                    {
                        for (long j = 1; j < blockheight - 1; j++) //for every row (except outside rows)
                        {
                            //Get the pixel values of center cell and  neighboring cells

                            center = inpixels.GetValue(i, j);

                            topleft = inpixels.GetValue(i - 1, j + 1);
                            topmiddle = inpixels.GetValue(i, j + 1);
                            topright = inpixels.GetValue(i + 1, j + 1);
                            middleleft = inpixels.GetValue(i - 1, j);
                            middleright = inpixels.GetValue(i + 1, j);
                            bottomleft = inpixels.GetValue(i - 1, j - 1);
                            bottommiddle = inpixels.GetValue(i, j - 1);
                            bottomright = inpixels.GetValue(i - 1, j - 1);


                            //compare center cell value with middle left cell and middle right cell in a 3x3 grid. If true, give output raster value of 1
                            if ((Convert.ToDouble(center) >= Convert.ToDouble(middleleft)) && (Convert.ToDouble(center) >= Convert.ToDouble(middleright)))
                            {
                                outpixels.SetValue(1, i, j);
                            }


                            //compare center cell value with top middle and bottom middle cell in a 3x3 grid. If true, give output raster value of 1
                            else if ((Convert.ToDouble(center) >= Convert.ToDouble(topmiddle)) && (Convert.ToDouble(center) >= Convert.ToDouble(bottommiddle)))
                            {
                                outpixels.SetValue(1, i, j);
                            }

                            //if neither conditions are true, give raster value of 0
                            else
                            {

                                outpixels.SetValue(0, i, j);
                            }
                        }
                    }
                    //Write the pixel array to the pixel block.
                    outpixelblock3.set_PixelData(k, outpixels);
                }
                //Finally, write the pixel block back to the raster.
                tlc = outRasterCursor.TopLeft;
                rasterEdit.Write(tlc, (IPixelBlock)outpixelblock3);
            }
            while (inrasterCursor.Next() == true && outRasterCursor.Next() == true);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);


        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message);
        }

    }
    
respondido por el Conor 26.09.2013 - 01:06
10

Veo que esto ya se resolvió con el Póster Original (OP), pero publicaré una solución simple en Python en caso de que alguien en el futuro esté interesado en diferentes formas de resolver este problema. Soy parcial al software de código abierto, así que aquí hay una solución que usa 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 ]]

Luego, itera a través de tus datos con un bucle anidado:

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

O tal vez desee aplanar su matriz 2-D con una lista de comprensión:

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

De todos modos, mientras se iteran a través de los datos celda por celda, es posible incluir algunos condicionales en su bucle para cambiar / editar valores. Vea este script que escribí para diferentes formas de acceder a los datos: enlace .

    
respondido por el asonnenschein 27.09.2013 - 01:04
8

Actualizar! La solución numpy:

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

Por lo tanto, hacer que la matriz terminada vuelva al ráster usando arcpy es problemático. arcpy.NumPyArrayToRaster es squirrelly y tiende a redefinir extensiones incluso si le das tus coordenadas LL.

Prefiero guardar como texto.

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

Estoy ejecutando Python como 64 bits para la velocidad; a partir de ahora, esto significa que no puedo enviar un encabezado a numpy.savetxt. Así que tengo que abrir la salida y agregar el encabezado ASCII que Arc quiere, antes de convertir ASCII a Raster

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

La versión numpy ejecuta mi ráster de turnos, multiplicaciones y suma mucho más rápido (1000 iteraciones en 2 minutos) que la versión arcpy (1000 iteraciones en 15 min)

VERSIÓN ANTIGUA Puedo eliminar esto más tarde Acabo de escribir un guión similar. Intenté convertir a puntos y usar el cursor de búsqueda. Obtuve solo 5000 iteraciones en 12 horas. Por lo tanto, busqué otra manera.

Mi forma de hacerlo es recorrer las coordenadas del centro de la celda de cada celda. Comienzo en la esquina superior izquierda y me muevo de derecha a izquierda. Al final de la fila, me muevo hacia abajo una fila y comienzo de nuevo a la izquierda. Tengo un ráster de 240 m con 2603 columnas y 2438 filas, por lo que hay un total de 6111844 celdas en total. Yo uso una variable iterador y un bucle while. Vea abajo

Algunas notas: 1 - necesitas saber las coordenadas de la extensión

2: ejecute con coordenadas de punto para el centro de la celda: mueva a 1/2 de tamaño de celda desde los valores de extensión

3 - Mi secuencia de comandos está usando el valor de celda para extraer un ráster específico del valor, luego cambiar este ráster al centro de la celda original. Esto se agrega a un ráster cero para expandir la extensión antes de agregarlo a un ráster final. Esto es solo un ejemplo. Puedes poner tus sentencias condicionales aquí (segunda sentencia if dentro del bucle while).

4 - Este script asume que todos los valores ráster se pueden convertir como enteros. Esto significa que necesita deshacerse de los datos sin datos primero. Con IsNull.

6 - Todavía no estoy contento con esto y estoy trabajando para sacarlo completamente de arcpy. Preferiría crear matrices numpy y hacer los cálculos allí y luego volver a Arc.

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")
    
respondido por el SamEPA 27.09.2013 - 23:50
4

Intente usar IGridTable, ICursor, IRow. Este fragmento de código es para actualizar valores de celdas ráster, sin embargo, muestra los conceptos básicos de iteración:

¿Cómo puedo agregar un nuevo campo en una tabla de atributos de ráster y recorrerlo?

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é atravesando la tabla, puede obtener el valor de la fila del campo específico usando row.get_Value(yourfieldIndex) . Si tu Google

  

arcobjects row.get_Value

debería poder obtener muchos ejemplos que muestran esto.

Espero que ayude.

    
respondido por el artwork21 23.09.2013 - 03:40
4

¿Qué tal esto como una idea radical? Requeriría que programes en Python o ArcObjects.

  1. Convierta su cuadrícula en una clase de puntos.
  2. Crea campos XY y rellena.
  3. Cargue los puntos en un diccionario donde key es una cadena de X, Y y item es el valor de la celda ..
  4. Recorra su diccionario y, para cada punto, resuelva los 8 XY de las celdas circundantes.
  5. Recupéralas de tu diccionario y prueba con tus reglas, tan pronto como encuentres un valor verdadero, puedes omitir el resto de las pruebas.
  6. Escriba los resultados en otro diccionario y luego vuelva a convertirlos en una cuadrícula creando primero una FeatureClass de puntos y luego convierta los puntos en una cuadrícula.
respondido por el Hornbydd 24.09.2013 - 18:49
1

Los datos ráster AFAIK se pueden leer de tres maneras:

  • por celda (ineficiente);
  • por imagen (bastante eficiente);
  • por bloques (la forma más eficiente).

Sin reinventar la rueda, sugiero leer estos diapositivas de Chris Garrard.

Por lo tanto, el método más eficiente es leer los datos por bloque, sin embargo, esto causaría una pérdida de datos en la correspondencia de los píxeles ubicados sobre los límites del bloque al aplicar el filtro. Por lo tanto, una forma alternativa segura debería consistir en leer toda la imagen a la vez y utilizar el método numpy.

En el lado computacional, en cambio, debería usar gdalfilter.py e implícitamente el enfoque VRT KernelFilteredSource para aplicar los filtros necesarios y, sobre todo, evitar cálculos pesados.

    
respondido por el Antonio Falciano 29.09.2013 - 10:29

Lea otras preguntas en las etiquetas