3

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.

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]]
  • Sorry, I found out a bug. I'm fixing it. – xunilk Feb 22 '17 at 12:18
  • I fixed the bug in my code. – xunilk Feb 22 '17 at 12:24
2

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.

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.