Python – Buffering Pixels in an Array

gdalnumpypythonscipy

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.

Best Answer

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]]
Related Question