Buffering pixels in an array- python

Question

Is there a way in SciPy or NumPy to buffer values in an array?

I have several rasters that I read as arrays using gdal to do some math/masking and then I write the final array back to geoTIFF. One of the layers I use to make my outputs needs to have the 0's in the array be buffered by one pixel before it can be used as a mask. Is there a way to do this to an array in Python? So if my input array looked like this:

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

Then the out array would look like this:

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

If there's not a straightforward way to do this with arrays, what's the best way to do this in GDAL? I only want the 0's in the raster/array to be buffered with 0's, and just by one pixel.

Answer 1

It's possible to buffer without 'for loops' with one line of code.

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. To buffer 0-value pixels:

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

2. To buffer 1-value pixels:

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

Answer 2

Next code works for your array example.

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)

After running it at the Python Console of QGIS I got your desired array:

[[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 3

You can define a simple function to return the slice of array values which are adjacent to a given index...

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]

Then you can use this function to modify a copy of your array, changing the value at each index if the value you are buffering is present in this adjacent slice:

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

It would probably be best to wrap this second part into a separate function too.

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

This method also allows you to easily buffer your raster with a different buffer distance, by changing the dist argument.

Answer 4

While the previous answers worked for a small array, when dealing with my 4800 x 4800 arrays, the solution was way too slow. The following function works in my case because I am buffering the 0-value pixels. It could be manipulated to work with other values, though:

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

This could also probably be broken down to even smaller functions as a lot of it is repetitive, but in the interest of time this works as well.