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

13

Vea este enlace para más detalles.

El problema:

Quiero recorrer un ráster continuo (uno que no tiene tabla de atributos), celda por celda, y obtener el valor de la celda. Quiero tomar esos valores y ejecutar condicionales sobre ellos, emulando los pasos de álgebra de mapas detallados a continuación sin usar la calculadora ráster.

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

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

El análisis necesitaba:

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

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

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

Donde el archivo 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, proporcione el valor de 1:

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

Donde el archivo 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 las celdas superior e inferior derecha, dé el valor de 1:

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

Donde el archivo 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 las celdas de abajo a la izquierda y hacia arriba, dé el valor de 1:

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

Donde el archivo kernel se ve así:

3 3 
0 0 1
0 0 0
1 0 0 

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

¿Por qué no usar álgebra de mapas?

A continuación se ha señalado que mi problema podría resolverse utilizando el álgebra de mapas, pero como se ve 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 individualmente siete veces y utilizar un poco de memoria para crear siete rásteres.

¿Cómo se debe atacar el problema?

El enlace anterior aconseja utilizar la interfaz IPixelBlock, sin embargo, de la documentación de ESRI no está claro si realmente está accediendo a un solo valor de celda a través de IPixelBlock, o si está accediendo a múltiples valores de celda desde el tamaño del IPixelBlock que configuró. 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 obvio.

En resumen:

¿Cuál es el mejor método para recorrer cada celda en un ráster CONTINUO (que no tiene una 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.

Conor
fuente
44
Casi siempre es innecesario recorrer cada celda de una trama. ¿Puede proporcionar más información sobre lo que está tratando de hacer?
usuario2856
2
@Luke es correcto: con mucho, la mejor manera de realizar un cálculo de ráster iterativo en cualquier SIG es evitar el bucle explícito a través de las celdas, porque debajo del capó ya se ha optimizado el bucle que se debe hacer. En su lugar, busque una manera de utilizar la funcionalidad de álgebra de mapas proporcionada por el SIG, si eso es posible. Si tuviera que describir su análisis, podría obtener respuestas útiles que utilicen dicho enfoque.
whuber
@Luke He agregado detalles del análisis.
Conor
1
Gracias por la aclaración, Conor. Estoy de acuerdo en que si su SIG incurre en una sobrecarga considerable para cada cálculo de trama, escribir su propio bucle puede ser más eficiente. Por curiosidad, ¿cuál es la interpretación prevista de este (inusual) conjunto de condiciones?
whuber
1
@whuber Es para operaciones de detección de bordes crear polígonos vectoriales a partir de mi ráster. La aplicación es conceptualmente similar a la identificación de cuencas hidrológicas de un DEM (piense en la celda central en las estadísticas del vecindario enumeradas anteriormente como el "pico" del que fluiría el agua aguas abajo) pero está fuera del campo de la hidrología. Anteriormente he estado usando Flow Direction y Basin Rasters para este propósito, pero son propensos a errores en mi análisis final debido a que las propiedades de estos métodos no son exactamente lo que necesito.
Conor

Respuestas:

11

Veo que esto ya ha sido resuelto por el Cartel 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 partidario del 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)

Implemente la función de esta manera:

#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, repita sus 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 comprensión de la lista:

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

De todos modos, mientras se itera a través de los datos celda por celda, es posible lanzar algunos condicionales en su ciclo para cambiar / editar valores. Consulte este script que escribí para conocer las diferentes formas de acceder a los datos: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py .

asonnenschein
fuente
Me gusta la simplicidad y elegancia de esta solución. Esperaré unos días más y, si nadie más llega con una solución de igual o mayor calidad, agregaré etiquetas para ampliar el alcance de la pregunta en beneficio de la comunidad y otorgarle la recompensa.
Conor
Gracias, @Conor! Nos encontramos con un problema similar en mi lugar de trabajo a principios de esta semana, así que lo resolví escribiendo una clase con GDAL / python. Específicamente, necesitábamos un método del lado del servidor para calcular el valor medio de un área de un ráster dado solo un cuadro delimitador de un usuario en nuestra aplicación del lado del cliente. ¿Crees que sería beneficioso si agregara el resto de la clase que escribí?
asonnenschein
Sería útil agregar código que muestre cómo leer la matriz 2D que recuperó y editar sus valores.
Conor
9

¡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 a la trama usando arcpy es problemático. arcpy.NumPyArrayToRaster es ardilla y tiende a redefinir extensiones incluso si lo alimenta a sus coordenadas LL.

Prefiero guardar como texto.

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

Estoy ejecutando Python como velocidad de 64 bits; en este momento, esto significa que no puedo alimentar numpy.savetxt un encabezado. Entonces tengo que abrir la salida y agregar el encabezado ASCII que Arc quiere antes de convertir ASCII a Ráster

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 turno, multiplicaciones y sumas mucho más rápido (1000 iteraciones en 2 minutos) que la versión arcpy (1000 iteraciones en 15 min)

VERSIÓN ANTERIOR Puedo eliminar esto más tarde, acabo de escribir un script similar. Intenté convertir a puntos y usar el cursor de búsqueda. Obtuve solo 5000 iteraciones en 12 horas. Entonces, busqué otra forma.

Mi forma de hacer esto es iterar a través de las coordenadas del centro celular 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 y empiezo de nuevo a la izquierda. Tengo un ráster de 240 m con 2603 columnas y 2438 filas, por lo que un total de 6111844 celdas totales. Yo uso una variable iteradora y un bucle while. Vea abajo

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

2: ejecutar con coordenadas de punto para el centro de la celda: mover 1/2 del tamaño de la celda desde los valores de extensión

3 - Mi script está usando el valor de celda para extraer un ráster específico de valor, luego desplazar este ráster para centrarlo en 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. Puede poner sus declaraciones condicionales aquí (segunda declaración if dentro del ciclo while).

4: este script supone que todos los valores ráster se pueden convertir como enteros. Esto significa que primero debe deshacerse de la ausencia de datos. Con IsNull.

6 - Todavía no estoy contento con esto y estoy trabajando para sacarlo de Arcpy por completo. Prefiero lanzar como matrices numpy y hacer los cálculos allí y luego llevarlo de vuelta 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")
revs SamEPA
fuente
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 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 fila de campo específico mediante el uso row.get_Value(yourfieldIndex). Si usted google

arcobjects row.get_Value

deberías poder obtener muchos ejemplos que lo demuestren.

Espero que ayude.

artwork21
fuente
1
Desafortunadamente, no tuve en cuenta, y editaré en mi pregunta original anterior, que mi ráster tiene muchos valores continuos que consisten en valores dobles grandes, y como tal, este método no funcionará porque mi ráster no tiene valores de tabla de atributos.
Conor
4

¿Qué tal esto como una idea radical? Te requeriría programar en Python o ArcObjects.

  1. Convierta su cuadrícula en un punto de clases de características.
  2. Crea campos XY y rellena.
  3. Cargue los puntos en un diccionario donde la clave es una cadena de X, Y y el elemento es el valor de la celda.
  4. Revise su diccionario y para cada punto calcule los 8 XY de la celda circundante.
  5. Recupere estos de su diccionario y pruebe con sus reglas, tan pronto como encuentre un valor que sea verdadero, puede omitir el resto de las pruebas.
  6. Escriba los resultados en otro diccionario y luego vuelva a convertirlos en una cuadrícula creando primero un FeatureClass de puntos y luego convierta los puntos en una cuadrícula.
Hornbydd
fuente
2
Al convertir a un conjunto de características de puntos, esta idea elimina las dos cualidades de la representación de datos basada en ráster que la hacen tan efectiva: (1) encontrar vecinos es una operación extremadamente simple en tiempo constante y (2) porque el almacenamiento explícito de ubicaciones es no necesarios, los requisitos de RAM, disco y E / S son mínimos. Por lo tanto, aunque este enfoque funcionará, es difícil encontrar alguna razón para recomendarlo.
whuber
Gracias por tu respuesta Hornbydd. Estoy de acuerdo con implementar un método como este, pero parece que los pasos 4 y 5 no serían muy efectivos desde el punto de vista computacional. Mis rásteres tendrán un mínimo de 62,500 celdas (la resolución mínima para mi ráster que he establecido es de 250 celdas x 250 celdas, pero la resolución puede y generalmente consiste en mucho más), y tendría que hacer una consulta espacial para cada condición para ejecutar mis comparaciones ... Como tengo 6 condiciones, eso sería 6 * 62500 = 375000 consultas espaciales. Estaría mejor con álgebra de mapas. Pero gracias por esta nueva forma de ver el problema. Votado
Conor
¿No puedes convertirlo a ASCII y luego usar un programa como R para hacer la computación?
Oliver Burdekin
Además, tengo un applet de Java que escribí que podría modificarse fácilmente para satisfacer sus condiciones anteriores. Era solo un algoritmo de suavizado, pero las actualizaciones serían bastante fáciles de hacer.
Oliver Burdekin
Siempre que se pueda llamar al programa desde la plataforma .NET para un usuario que solo tenga instalado .NET Framework 3.5 y ArcGIS 10. El programa es de código abierto y pretendo que esos sean los únicos requisitos de software cuando se entregan a los usuarios finales. Si su respuesta se puede implementar para cumplir estos dos requisitos, se considerará una respuesta válida. Agregaré una etiqueta de versión a la pregunta también para aclaración.
Conor
2

Una solución:

Resolví esto más temprano 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 para interactuar con la trama. El siguiente método 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 hace las 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 en el ráster. Cuando crea IPixelBlock, le indica 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, configúrelo como sus parámetros IPixelBlock. Si desea recorrer todo el ráster, debe configurar IPixelBlock igual al tamaño de todo el ráster. Hago esto en el siguiente método pasando el tamaño a IRasterCursor (pSize) y luego obteniendo el PixelBlock del cursor ráster.

La otra clave es que debe 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, escriba todo su SafeArray de nuevo en IPixelBlock, luego escriba su IPixelBlock en IRasterCursor, y finalmente use IRasterCursor para establecer la ubicación para comenzar a escribir e IRasterEdit para hacer la escritura. Este último paso es donde realmente edita 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);
        }

    }
Conor
fuente
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 estas diapositivas esclarecedoras de Chris Garrard.

Entonces, 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 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 la imagen completa de una vez y usar el enfoque 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.

Antonio Falciano
fuente