Buffering pixels in an array- python

Question

¿Existe alguna forma en SciPy o NumPy de almacenar valores en un array?

Tengo varios rasters que leo como matrices utilizando gdal para hacer algunas matemáticas / enmascaramiento y luego escribo la matriz final de nuevo a geoTIFF. Una de las capas que utilizo para hacer mis salidas tiene que tener el 0 en la matriz de ser amortiguado por un píxel antes de que pueda ser utilizado como una máscara. ¿Hay alguna manera de hacer esto a una matriz en Python? Así que si mi matriz de entrada se parecía a esto:

1 0 1 1 1 1
0 0 1 1 1 1
1 0 1 1 1 1
1 1 1 0 1 1
1 1 1 1 1 1

Entonces la matriz de salida tendría este aspecto:

0 0 0 1 1 1
0 0 0 1 1 1
0 0 0 0 0 1
0 0 0 0 0 1
1 1 0 0 0 1

Si no hay una forma sencilla de hacerlo con arrays, ¿cuál es la mejor manera de hacerlo en GDAL? Solo quiero que los 0's en el raster/array sean bufferizados con 0's, y solo por un pixel.

Answer 1

El siguiente código funciona para tu ejemplo de array.

import numpy as np

#array is a matrix but afterward is easy to convert it in array with numpy

array = [[1, 0, 1, 1, 1, 1],
         [0, 0, 1, 1, 1, 1],
         [1, 0, 1, 1, 1, 1],
         [1, 1, 1, 0, 1, 1],
         [1, 1, 1, 1, 1, 1]]

indexes = []

for i, item in enumerate(array):
    for j, element in enumerate(item):
        if element == 0:
            if i-1 >= 0 and j-1 >= 0:
                tuple = (i-1, j-1)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i-1 >= 0 and j >= 0:
                tuple = (i-1, j)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i-1 >= 0 and j+1 >= 0:
                tuple = (i-1, j+1)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i >= 0 and j-1 >= 0:
                tuple = (i, j-1)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i >= 0 and j >= 0:
                tuple = (i, j)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i >= 0 and j+1 >= 0:
                tuple = (i, j+1)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i+1 >= 0 and j-1 >= 0:
                tuple = (i+1, j-1)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i+1 > 0 and j > 0:
                tuple = (i+1, j)
                if tuple not in indexes:
                    indexes.append(tuple)

            if i+1 > 0 and j+1 > 0:
                tuple = (i+1, j+1)
                if tuple not in indexes:
                    indexes.append(tuple)

buffered_array = [[1 for i in range(len(array[0]))] for i in range(len(array))]

for index in indexes:
    buffered_array[index[0]][index[1]] = 0

print np.array(buffered_array)

Después de ejecutarlo en la consola Python de QGIS obtuve el array deseado:

[[0 0 0 1 1 1]
 [0 0 0 1 1 1]
 [0 0 0 0 0 1]
 [0 0 0 0 0 1]
 [1 1 0 0 0 1]]

Answer 2

Puedes definir una función sencilla para devolver la porción de valores de la matriz que son adyacentes a un índice dado...

import numpy as np

def adjacent_slice(arr, index, dist=1):

    '''array is 2d numpy array, index is (row,col) tuple, dist is the 
       maximum buffer distance'''

    # if either index value - dist is less than 0, 0 should be used as
    # the start index
    row_start = max(0, index[0] - dist)
    col_start = max(0, index[1] - dist)

    row_stop = index[0] + dist + 1
    col_stop = index[1] + dist + 1

    return arr[row_start:row_stop, col_start:col_stop]

A continuación, puede utilizar esta función para modificar una copia de su matriz, cambiando el valor en cada índice si el valor que está almacenando está presente en esta rebanada adyacente:

arr = np.array([[1, 0, 1, 1, 1, 1],
                [0, 0, 1, 1, 1, 1],
                [1, 0, 1, 1, 1, 1],
                [1, 1, 1, 0, 1, 1],
                [1, 1, 1, 1, 1, 1]])

# do not modify in place!
output_arr = np.copy(arr)

for i in range(arr.shape[0]):
    for j in range(arr.shape[1]):
        if 0 in adjacent_slice(arr, (i,j)):
            output_arr[i,j] = 0

Probablemente también sería mejor incluir esta segunda parte en una función independiente.

def buffer_value(arr, buffer_val, dist=1):

    output_array = np.copy(arr)

    for i in range(arr.shape[0]):
        for j in range(arr.shape[1]):
            if buffer_val in adjacent_slice(arr, (i,j), dist):
                output_arr[i,j] = buffer_val

    return output_array

Este método también le permite almacenar fácilmente su trama con una distancia de almacenamiento diferente, cambiando el argumento dist.

Answer 3

Mientras que las respuestas anteriores funcionaban para una matriz pequeña, cuando se trataba de mis matrices de 4800 x 4800, la solución era demasiado lenta. La siguiente función funciona en mi caso porque estoy almacenando en búfer los píxeles de valor 0. Podría ser manipulado para trabajar con otros valores, sin embargo:

def bufferArray(Arr):
    nrows= Arr.shape[0] # number of rows (size of Y direction)
    ncols= Arr.shape[1] # number of columns (size of X)

    ones_row = np.ones((1, ncols), dtype = np.int) # add to top and bottom
    ones_col = np.ones((nrows), dtype = np.int) # add to left and right

    # Arr1 - shift up (remove first row, add bottom row)
    Arr1 = np.delete(Arr, (0), axis=0) # removes first row
    Arr1 = np.insert(Arr1, (nrows-1), ones_row, axis=0) # zero adds to the row axis

    # Arr2 - shift down (remove last row, add first row)
    Arr2 = np.delete(Arr, (nrows-1), axis=0) # removes last row
    Arr2 = np.insert(Arr2, 0, ones_row, axis=0)

    # Arr3 - shift right (remove last column, add first column)
    Arr3 = np.delete(Arr, (ncols-1), axis=1)
    Arr3 = np.insert(Arr3, 0, ones_col, axis=1)

    #Arr4 - shift left (remove first column, add last column)
    Arr4 = np.delete(Arr, (0), axis=1)
    Arr4 = np.insert(Arr4, (ncols-1), ones_col, axis=1)

    # Arr5 shift up one AND left one, Arr1 is already shifted up, now shift left
    Arr5 = np.delete(Arr1, (0), axis=1)
    Arr5 = np.insert(Arr5, (ncols-1), ones_col, axis=1)

    # Arr6 shift up one AND right one -- shift Arr1 right
    Arr6 = np.delete(Arr1, (ncols-1), axis=1)
    Arr6 = np.insert(Arr6, 0, ones_col, axis=1)

    # Arr7 shift down one AND left one -- shift Arr2 left
    Arr7 = np.delete(Arr2, (0), axis=1)
    Arr7 = np.insert(Arr7, (ncols-1), ones_col, axis=1)

    # Arr8 shift down one AND right one -- shift Arr2 right
    Arr8 = np.delete(Arr2, (ncols-1), axis=1)
    Arr8 = np.insert(Arr8, 0, ones_col, axis=1)

    bufferedArr = Arr*Arr1*Arr2*Arr3*Arr4*Arr5*Arr6*Arr7*Arr8

    return bufferedArr

Probablemente también podría dividirse en funciones aún más pequeñas, ya que muchas son repetitivas, pero en aras del tiempo esto también sirve.

Answer 4

Es posible almacenar en búfer sin 'bucles for' con una línea de código.

import numpy as np
arr = np.array([[1, 0, 1, 1, 1, 1],
                [0, 0, 1, 1, 1, 1],
                [1, 0, 1, 1, 1, 1],
                [1, 1, 1, 0, 1, 1],
                [1, 1, 1, 1, 1, 1]])

1. Para almacenar píxeles de valor 0:

   from scipy.ndimage import minimum_filter
   buffer_size = 1
   minimum_filter(arr, size=2*buffer_size+1, mode='constant', cval=1)

2. Para almacenar píxeles de 1 valor:

   from scipy.ndimage import maximum_filter
   buffer_size = 1
   maximum_filter(arr, size=2*buffer_size+1, mode='constant', cval=0)