Combine... data classification from raster combinations

[DanPatterson_Retired]

I rarely get a great idea.  This is another one of those cases.  But I do recognize a great idea when I see it.  Sadly I usually forget it, then rediscover it many months or years later.  I spotted a fantastic idea today and I had to drop everything (between answering questions) and jot it down.

The premise that twigged my interest was about how to 'combine' values of elevation, slope and aspect on Stack Overflow.  I did a wowwser!  Terrain derivatives on Stack!  I was enthralled.  The author.... unutbu ... of mysterious persona, put together a short missive using built-in numpy tools.  The code is short and sweet, works with floating point or integer rasters (aka, arrays) and is fast (at least in my testing so far).

Before I ruin the rest of my day with some droll corner case I hadn't thought of, I will jot this down now for posterity.

The large print

• This doesn't require a Spatial Analyst License
• A complete implementation (coming soon) uses arcpy's RasterToNumPyArray and NumPyArrayToRaster to communicate between the thin veil separating Arc* from the array world
• Rasters/arrays with nodata values will be implemented
• This is in the initial stages of development of a full-fledged tool, so hold you ...but!... comments

Start with 3 small rasters/arrays.  This way, you can see how the system works.

Array 'a'	Array 'b'
array([[0, 0, 0, 4, 4, 4, 1, 1, 1],        [0, 0, 0, 4, 4, 4, 1, 1, 1],        [0, 0, 0, 4, 4, 4, 1, 1, 1],        [2, 2, 2, 1, 1, 1, 2, 2, 2],        [2, 2, 2, 1, 1, 1, 2, 2, 2],        [2, 2, 2, 1, 1, 1, 2, 2, 2],        [1, 1, 1, 4, 4, 4, 0, 0, 0],        [1, 1, 1, 4, 4, 4, 0, 0, 0],        [1, 1, 1, 4, 4, 4, 0, 0, 0]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍	array([[0, 0, 0, 1, 1, 1, 2, 2, 2],        [0, 0, 0, 1, 1, 1, 2, 2, 2],        [0, 0, 0, 1, 1, 1, 2, 2, 2],        [3, 3, 3, 4, 4, 4, 5, 5, 5],        [3, 3, 3, 4, 4, 4, 5, 5, 5],        [3, 3, 3, 4, 4, 4, 5, 5, 5],        [0, 0, 0, 1, 1, 1, 2, 2, 2],        [0, 0, 0, 1, 1, 1, 2, 2, 2],        [0, 0, 0, 1, 1, 1, 2, 2, 2]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Array 'c'	unique combinations
array([[0, 0, 0, 3, 3, 3, 0, 0, 0],        [0, 0, 0, 3, 3, 3, 0, 0, 0],        [0, 0, 0, 3, 3, 3, 0, 0, 0],        [1, 1, 1, 4, 4, 4, 1, 1, 1],        [1, 1, 1, 4, 4, 4, 1, 1, 1],        [1, 1, 1, 4, 4, 4, 1, 1, 1],        [2, 2, 2, 5, 5, 5, 2, 2, 2],        [2, 2, 2, 5, 5, 5, 2, 2, 2],        [2, 2, 2, 5, 5, 5, 2, 2, 2]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍	array([[0, 0, 0, 6, 6, 6, 1, 1, 1],        [0, 0, 0, 6, 6, 6, 1, 1, 1],        [0, 0, 0, 6, 6, 6, 1, 1, 1],        [2, 2, 2, 7, 7, 7, 3, 3, 3],        [2, 2, 2, 7, 7, 7, 3, 3, 3],        [2, 2, 2, 7, 7, 7, 3, 3, 3],        [4, 4, 4, 8, 8, 8, 5, 5, 5],        [4, 4, 4, 8, 8, 8, 5, 5, 5],        [4, 4, 4, 8, 8, 8, 5, 5, 5]],‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So here it is.

def combine_(*arrs):
    """Combine arrays to produce a unique classification scheme
    :
    :References:
    :-----------
    : https://stackoverflow.com/questions/48035246/
    :       intersect-multiple-2d-np-arrays-for-determining-zones
    : original 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, uniqs‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Now stop reading here if you aren't interested in the details, but have a quasi-firm grasp on the process.  Read on if you must, but suffice to say unubtu was on to something.

-----------------------------  details .... details .... details... read at your peril  ------------------------------------------------------------

In essence, what the code does, is find the unique values of each array shoved into the function.  In our case, the list of unique values is short (see line 3 below).  Not too onerous a combination slate.

Line 10 can be modified to show the unique values in each array... useful if they aren't sequential.  This is show below in line 1.

induni = [np.unique(arr, return_inverse=True)[0] for arr in arrs]  # arrays a, b, c

[array([0, 1, 2, 4]), array([0, 1, 2, 3, 4, 5]), array([0, 1, 2, 3, 4, 5])]‍‍‍‍‍‍‍‍‍‍‍‍

The nifty bit, is how unutbu constructs the 'indices' in two steps, using the maximum value in each array, determines the unique combination, applies it to the arrays to get a unique value ('val') for each combination and then reassigns a class ('label' to it).

This process, not documented in the code, is expanded upon here.  In lines 1-4 below, an array, combos, is constructed from the values in the array plus the 'val' derived from the above process.

idx = [np.unique(arr, return_inverse=True)[1] for arr in arrs]
vals = (indices * strides).sum(axis=-1)
dt = [('a', '<i8'), ('b', '<i8'), ('c', '<i8'), ('vals', '<i8')]
combos = np.array(list(zip(*idx)), dtype= dt)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

If the 'combos' array is sorted, the unique values can be extracted from it.  The values in the 'vals' column is used to assign a 'label' to each unique combination.

srt_order = np.argsort(combos, order=['a','b', 'c'])  # sort the combos

srted = combos[srt_order]  # ---- extract them in sorted order

np.unique(srted)  # ---- determine the unique combinations
 
array([(0, 0, 0,   0),
       (0, 2, 2,  56),
       (1, 0, 2,  49),
       (1, 2, 0,   9),
       (1, 4, 4, 113),
       (2, 3, 1,  38),
       (2, 5, 1,  46),
       (3, 1, 3,  79),
       (3, 1, 5, 127)],
      dtype=[('a', '<i8'), ('b', '<i8'),
             ('c', '<i8'), ('vals', '<i8')])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Hope it is clear as mud, and I have some experimenting to do, but I just thought I would throw it out there.