Intersect multiple 2D np arrays for determining zones

Question

Using this small reproducible example, I've so far been unable to generate a new integer array from 3 arrays that contains unique groupings across all three input arrays.

The arrays are related to topographic properties:

import numpy as np
asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

The idea is that the geographic contours are broken into 3 different properties using GIS routines:

1-8 for aspect (1=north facing, 2=northeast facing, etc.)
9-12 for slope (9=gentle slope...12=steepest slope)
13-16 for elevation (13=lowest elevations...16=highest elevations)

The small graphic below attempts to depict the kind of result I'm after (array shown in lower left). Note, the "answer" given in the graphic is but one possible answer. I'm not concerned about the final arrangement of integers in the resulting array so long as the final array contains an integer at each row/column index that identifies unique groupings.

For example, the array indexes at [0,1] and [0,2] have the same aspect, slope, and elevation and therefore receive the same integer identifier in the resulting array.

Does numpy have a built in routine for this kind of thing?

Although Numpy doesn't have a groupby functioon, there's a library which provides an efficient one, according to its author. See https://stackoverflow.com/a/38015063/4014959 — PM 2Ring, Dec 30 '17 at 16:44
See a related follow-up post about spatially contiguous zones here: https://stackoverflow.com/questions/54602530/intersect-multiple-2d-np-arrays-for-determining-spatially-contiguous-zones — user2256085, Feb 09 '19 at 05:04

unutbu · Accepted Answer · 2017-12-30T19:50:11.707

Each location in the grid is associated with a tuple composed of one value from asp, slp and elv. For example, the upper left corner has tuple (8,9,13). We would like to map this tuple to a number which uniquely identifies this tuple.

One way to do that would be to think of (8,9,13) as the index into the 3D array np.arange(9*13*17).reshape(9,13,17). This particular array was chosen to accommodate the largest values in asp, slp and elv:

In [107]: asp.max()+1
Out[107]: 9

In [108]: slp.max()+1
Out[108]: 13

In [110]: elv.max()+1
Out[110]: 17

Now we can map the tuple (8,9,13) to the number 1934:

In [113]: x = np.arange(9*13*17).reshape(9,13,17)

In [114]: x[8,9,13]
Out[114]: 1934

If we do this for each location in the grid, then we get a unique number for each location. We could end right here, letting these unique numbers serve as labels.

Or, we can generate smaller integer labels (starting at 0 and increasing by 1) by using np.unique with return_inverse=True:

uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

So, for example,

import numpy as np

asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

x = np.arange(9*13*17).reshape(9,13,17)
vals = x[asp, slp, elv]
uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

yields

array([[11,  0,  0,  1],
       [ 9, 12,  2,  3],
       [10,  8,  5,  3],
       [ 7,  6,  6,  4]])

The above method works fine as long as the values in asp, slp and elv are small integers. If the integers were too large, the product of their maximums could overflow the maximum allowable value one can pass to np.arange. Moreover, generating such a large array would be inefficient. If the values were floats, then they could not be interpreted as indices into the 3D array x.

So to address these problems, use np.unique to convert the values in asp, slp and elv to unique integer labels first:

indices = [ np.unique(arr, return_inverse=True)[1].reshape(arr.shape) for arr in [asp, slp, elv] ]
M = np.array([item.max()+1 for item in indices])
x = np.arange(M.prod()).reshape(M)
vals = x[indices]
uniqs, labels = np.unique(vals, return_inverse=True)
labels = labels.reshape(vals.shape)

which yields the same result as shown above, but works even if asp, slp, elv were floats and/or large integers.

Finally, we can avoid the generation of np.arange:

x = np.arange(M.prod()).reshape(M)
vals = x[indices]

by computing vals as a product of indices and strides:

M = np.r_[1, M[:-1]]
strides = M.cumprod()
indices = np.stack(indices, axis=-1)
vals = (indices * strides).sum(axis=-1)

So putting it all together:

import numpy as np

asp = np.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = np.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = np.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

def find_labels(*arrs):
    indices = [np.unique(arr, return_inverse=True)[1] for arr in arrs]
    M = np.array([item.max()+1 for item in indices])
    M = np.r_[1, M[:-1]]
    strides = M.cumprod()
    indices = np.stack(indices, axis=-1)
    vals = (indices * strides).sum(axis=-1)
    uniqs, labels = np.unique(vals, return_inverse=True)
    labels = labels.reshape(arrs[0].shape)
    return labels

print(find_labels(asp, slp, elv))

# [[ 3  7  7  0]
#  [ 6 10 12  4]
#  [ 8  9 11  4]
#  [ 2  5  5  1]]

a nice implementation and very speedy. It also makes histogram generation much easier. — NaN, Jan 19 '18 at 02:33

score 2 · Answer 2 · answered Dec 30 '17 at 16:08

This can be done using numpy.unique() and then a mapping like:

Code:

combined = 10000 * asp + 100 * slp + elv
unique = dict(((v, i + 1) for i, v in enumerate(np.unique(combined))))
combined_unique = np.vectorize(unique.get)(combined)

Test Code:

import numpy as np

asp = np.array([8, 1, 1, 2, 7, 8, 2, 3, 7, 6, 4, 3, 6, 5, 5, 4]).reshape((4, 4))  # aspect
slp = np.array([9, 10, 10, 9, 9, 12, 12, 9, 10, 11, 11, 9, 9, 9, 9, 9]).reshape((4, 4))  # slope
elv = np.array([13, 14, 14, 13, 14, 15, 16, 14, 14, 15, 16, 14, 13, 14, 14, 13]).reshape((4, 4))

combined = 10000 * asp + 100 * slp + elv
unique = dict(((v, i + 1) for i, v in enumerate(np.unique(combined))))
combined_unique = np.vectorize(unique.get)(combined)

print(combined_unique)

Results:

[[12  1  1  2]
 [10 13  3  4]
 [11  9  6  4]
 [ 8  7  7  5]]

asheets · Answer 3 · 2017-12-30T16:23:22.707

This seems like a similar problem to labeling unique regions in an image. This is a function I've written to do this, though you would first need to concatenate your 3 arrays to 1 3D array.

def labelPix(pix):
    height, width, _ = pix.shape
    pixRows = numpy.reshape(pix, (height * width, 3))
    unique, counts = numpy.unique(pixRows, return_counts = True, axis = 0)

    unique = [list(elem) for elem in unique]

    labeledPix = numpy.zeros((height, width), dtype = int)
    offset = 0
    for index, zoneArray in enumerate(unique):
        index += offset
        zone = list(zoneArray)
        zoneArea = (pix == zone).all(-1)
        elementsArray, numElements = scipy.ndimage.label(zoneArea)

        elementsArray[elementsArray!=0] += offset

        labeledPix[elementsArray!=0] = elementsArray[elementsArray!=0]

        offset += numElements

    return labeledPix

This will label unique 3-value combinations, while also assigning separate labels to zones which have the same 3-value combination, but are not in contact with one another.

asp = numpy.array([8,1,1,2,7,8,2,3,7,6,4,3,6,5,5,4]).reshape((4,4))  #aspect
slp = numpy.array([9,10,10,9,9,12,12,9,10,11,11,9,9,9,9,9]).reshape((4,4))  #slope
elv = numpy.array([13,14,14,13,14,15,16,14,14,15,16,14,13,14,14,13]).reshape((4,4)) #elevation

pix = numpy.zeros((4,4,3))
pix[:,:,0] = asp
pix[:,:,1] = slp
pix[:,:,2] = elv

print(labelPix(pix))

returns:

[[ 0  1  1  2]
 [10 12  3  4]
 [11  9  6  4]
 [ 8  7  7  5]]

PM 2Ring · Answer 4 · 2017-12-30T16:24:04.827

Here's a plain Python technique using itertools.groupby. It requires the input to be 1D lists, but that shouldn't be a major issue. The strategy is to zip the lists together, along with an index number, then sort the resulting columns. We then group identical columns together, ignoring the index number when comparing columns. Then we gather the index numbers from each group, and use them to build the final output list.

from itertools import groupby

def show(label, seq):
    print(label, ' '.join(['{:2}'.format(u) for u in seq]))

asp = [8, 1, 1, 2, 7, 8, 2, 3, 7, 6, 4, 3, 6, 5, 5, 4] 
slp = [9, 10, 10, 9, 9, 12, 12, 9, 10, 11, 11, 9, 9, 9, 9, 9] 
elv = [13, 14, 14, 13, 14, 15, 16, 14, 14, 15, 16, 14, 13, 14, 14, 13]

size = len(asp)
a = sorted(zip(asp, slp, elv, range(size)))
groups = sorted([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))
final = [0] * size
for i, g in enumerate(groups, 1):
    for j in g:
        final[j] = i

show('asp', asp)
show('slp', slp)
show('elv', elv)
show('out', final)

output

asp  8  1  1  2  7  8  2  3  7  6  4  3  6  5  5  4
slp  9 10 10  9  9 12 12  9 10 11 11  9  9  9  9  9
elv 13 14 14 13 14 15 16 14 14 15 16 14 13 14 14 13
out  1  2  2  3  4  5  6  7  8  9 10  7 11 12 12 13

There's no need to do that second sort, we could just use a plain list comp

groups = [[u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1])]

or generator expression

groups = ([u[-1] for u in g] for _, g in groupby(a, key=lambda t:t[:-1]))

I only did it so that my output matches the output in the question.

score 1 · Answer 5 · answered Dec 30 '17 at 17:17

Here's one way to solve this problem using a dictionary based lookup.

from collections import defaultdict
import itertools

group_dict = defaultdict(list)
idx_count = 0

for a, s, e in np.nditer((asp, slp, elv)):
     asp_tuple = (a.tolist(), s.tolist(), e.tolist())
     if asp_tuple not in group_dict:
         group_dict[asp_tuple] = [idx_count+1]
         idx_count += 1
     else:
         group_dict[asp_tuple].append(group_dict[asp_tuple][-1])


list1d = list(itertools.chain(*list(group_dict.values())))

np.array(list1d).reshape(4, 4)

# result 
array([[ 1,  2,  2,  3],
       [ 4,  5,  6,  7],
       [ 7,  8,  9, 10],
       [11, 12, 12, 13]])

Intersect multiple 2D np arrays for determining zones

5 Answers5

Code:

Test Code:

Results:

Linked