CODEWITHSUNDEEP
pythonnumpy
Jen
Jen

Reputation: 21

Finding the lowest neighboring index in a Python array

I have a problem in which we're supposed to write a function that, when given an input of a 2D array, will return the offset in rows and in columns of each index's neighboring index of the lowest value; one array for the offset of each index in rows and one array for the offset in columns. For example, if an index's lowest neighboring cell is down one row and to the right one column, the offset is 1,1; if the lowest neighboring cell is to the left, the offset is 0,-1; if it is the lowest cell out of its neighboring cells, then the offset is 0,0.

Because I couldn't find a faster and correct way to do this, I wrote a while loop that would iterate through each index and see which of the surrounding indexes for a point [i,j] was lower than all of the other surrounding indexes using a.all():

def findLowNhbr( terrain ):
    """Creates two 2D-arrays the shape of terrain consisting
    of the offsets (row and column) to the neighbor with the minimum eleveation"""
    rowOffset = np.zeros_like(terrain)
    colOffset = np.zeros_like(terrain)

for i in range(len(terrain)):
    if i == 0:
        rowOffset[i] = 0
        colOffset[i] = 0
    elif i == (len(terrain)-1):
        rowOffset[i] = 0
        colOffset[i] = 0
    else:
        for j in range(len(terrain[i])):
            if j == 0 or j == len(terrain[i])-1:
                rowOffset[:,j] = 0
                colOffset[:,j] = 0
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i-1,j-1]).all():
                rowOffset[i,j] = -1
                colOffset[i,j] = -1
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i,j-1]).all():
                rowOffset[i,j] = 0
                colOffset[i,j] = -1
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i+1,j-1]).all():
                rowOffset[i,j] = 1
                colOffset[i,j] = -1
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i-1,j]).all():
                rowOffset[i,j] = -1
                colOffset[i,j] = 0
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i+1,j]).all():
                rowOffset[i,j] = 1
                colOffset[i,j] = 0
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i-1,j+1]).all():
                rowOffset[i,j] = -1
                colOffset[i,j] = 1
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i,j]).all():
                rowOffset[i,j] = 0
                colOffset[i,j] = 1
            elif (terrain[i-1:i+2,j-1:j+2]>=terrain[i+1,j+1]).all():
                rowOffset[i,j] = 1
                colOffset[i,j] = 1
            else:
                rowOffset[i,j] = 0
                colOffset[i,j] = 0
return rowOffset, colOffset

It takes a long time to run, but it does run. I can't imagine that I'm actually doing this the most efficient way possible; any input?

Upvotes: 2

Views: 343

Answers (2)

Jaime
Jaime

Reputation: 67427

I like Mr. E's basic idea of stacking all the surrounding values in a single dimension, but I think there are better ways of creating the stacked array and converting the return of np.argmin to pairs of indices:

from numpy.lib.stride_tricks import as_strided

rows, cols = 100, 100
win_rows, win_cols = 3, 3 # these two should be odd numbers
terrain = np.random.rand(rows, cols)

# This takes a windowed view of the original array, no data copied
win_terrain = as_strided(terrain, shape=(rows-win_rows+1, cols-win_cols+1,
                                         win_rows, win_cols),
                         strides=terrain.strides*2)
# This makes a copy of the stacked array that will take up x9 times more memory
# than the original one
win_terrain = win_terrain.reshape(win_terrain.shape[:2] + (-1,))

indices = np.argmax(win_terrain, axis=-1)
offset_rows, offset_cols = np.unravel_index(indices,
                                            dims=(win_rows, win_cols))
# For some odd reason these arrays are non-writeable, so -= won't work
offset_rows = offset_rows - win_rows//2
offset_cols = offset_cols - win_cols//2

The resulting arrays are only (98, 98), i.e. the first and last column and rows are missing, since there isn't a fully defined window around them.

Upvotes: 1

YXD
YXD

Reputation: 32511

This should more or less do it in a vectorized manner, ignoring some issues at the boundaries for the moment which you can avoid by padding the input array with repeated values around the edges and trimming the output

import numpy as np

np.random.seed(0)
terrain = np.random.rand(10,10)

offsets = [(i,j) for i in range(-1,2) for j in range(-1,2)]

stacked = np.dstack( np.roll(np.roll(terrain,i,axis=0),j,axis=1) for i, j in offsets)

offset_index = np.argmin(stacked,axis=2)
output = np.array(offsets)[offset_index]

Explanation

  • Stack all the offsets into an NxMx9 array
  • Find the index of the minimum element (argmin) along this last axis axis=2
  • We turn this index into an array of offset vectors by using the result to index the offsets in the last line

Another perhaps cleaner way to get the all the initial offsets is:

from itertools import product
offsets = list(product((-1, 0, 1), (-1, 0, 1)))

Upvotes: 2

Related Questions