Sliding/Moving window operations in rasters and arrays

DanPatterson_Retired · ‎07-07-2016

Sliding/Moving windows ...

This is the companion to block functions introduced earlier. I will keep it simple. The stats functions for rasters with and without nodata values still apply to this type of treatment. The only difference is how the sub-arrays are generated. I won't get into the details of how they are done here, but the stride function uses some pretty tricky handling of the input array structure to parse out the moving windows.

Remember, block functions move one window size across and down, per movement cycle. Sliding functions move one cell across and down and sample the cell based on the window size...hence, there are many more windows to sample in sliding function than in block functions. I will just show the code snippets and you can play on your own. I have posted sample timing results for comparative purposes for both implementations.

Code section and comments

The stride_tricks does the work. You can examine the code in the ...site-packages, numpy lib folder.
functools provides decorator capabilities. I tend to import it by default since I use decorators a lot for timing and documentation functions.
The textwrap module provide several functions, with dedent being useful for output documentation.
The numpy set_printoptions function allows you to control the appearance of output. The version here is quite simplified...refer to the documentation for more information.

Stride and block functions provide the means to sample the input array according to the type of data, the size and layout. You must read the documentation and experiment in order to fully understand. It really messes with your head and makes the Rubic's cube look simple. The decorator and timer functions I have talked about before. Refer to my previous blogs for more information.

imports

import numpy as np
from numpy.lib.stride_tricks import as_strided
from functools import wraps
from textwrap import dedent
np.set_printoptions(edgeitems=3,linewidth=80,precision=2,
                    suppress=True,threshold=100)
‍‍‍‍‍‍

slide (stride) and block functions

#---- functions ----
def _check(a, r_c):
    """Performs the array checks necessary for stride and block.
    : a   - Array or list.
    : r_c - tuple/list/array of rows x cols.  
    :Attempts will be made to 
    :     produce a shape at least (1*c).  For a scalar, the
    :     minimum shape will be (1*r) for 1D array or (1*c) for 2D
    :     array if r<c.  Be aware
    """
    if isinstance(r_c, (int, float)):
        r_c = (1, int(r_c))
    r, c = r_c
    a = np.atleast_2d(a)
    shp = a.shape
    r, c = r_c = ( min(r, a.shape[0]), min(c, shp[1]) ) 
    a = np.ascontiguousarray(a)
    return a, shp, r, c, tuple(r_c)
    
def stride(a, r_c=(3, 3)):
    """Provide a 2D sliding/moving view of an array.  
    :  There is no edge correction for outputs.
    :
    :Requires
    :--------
    : a - array or list, usually a 2D array.  Assumes rows is >=1,
    :     it is corrected as is the number of columns.
    : r_c - tuple/list/array of rows x cols.  Attempts  to 
    :     produce a shape at least (1*c).  For a scalar, the
    :     minimum shape will be (1*r) for 1D array or 2D
    :     array if r<c.  Be aware
    """
    a, shp, r, c, r_c = _check(a, r_c)
    shape = (a.shape[0] - r + 1, a.shape[1] - c + 1) + r_c
    strides = a.strides * 2
    a_s = (as_strided(a, shape=shape, strides=strides)).squeeze()    
    return a_s


def block(a, r_c=(3, 3)):
    """See _check and/or stride for documentation.  This function
    :  moves in increments of the block size, rather than sliding
    :  by one row and column
    :
    """
    a, shp, r, c, r_c = _check(a, r_c)
    shape = (a.shape[0]/r, a.shape[1]/c) + r_c
    strides = (r*a.strides[0], c*a.strides[1]) + a.strides
    a_b = as_strided(a, shape=shape, strides=strides).squeeze()
    return a_b
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

timing decorator, decorated stat function and time test

def delta_time(func):
    """simple timing decorator function"""
    import time
    @wraps(func)
    def wrapper(*args, **kwargs):
        print("Timing function for... {}".format(func.__name__))
        t0 = time.perf_counter()        # start time
        result = func(*args, **kwargs)  # ... run the function ...
        t1 = time.perf_counter()        # end time
        print("Results for... {}".format(func.__name__))
        print("  time taken ...{:12.9e} sec.".format(t1-t0))
        #print("\n {}".format(result))  # print within wrapper
        return result                   # return result
    return wrapperdelta_time
def array_mean(a):
    """change the func"""
    a_mean = a.mean(axis=(2, 3))
    return None
def time_test():
    """time test for block and sliding raster"""
    N = [100, 500, 1000, 2000, 3000, 4000, 5000]
    for n in N:
        a = np.arange(n*n).reshape((n, n))
        #a0 = stride(a, block=(3, 3))  # uncomment this or below
        a0 = block(a, block=(3, 3))
        print("\nArray size... {0}x{0}".format(n))
        a_stat = array_mean(a0)       # time function
    return None
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

main section

if __name__=="__main__":
    """comparison between block and slide view"""
    prn=False
    time_test()
‍‍‍‍

Description: This code does...timing examples for sliding and block mean

Timing results

sliding mean block mean

array size total time total time (seconds)

1000x1000 3.85e-02 4.65...e-03

2000x2000 1.51e-01 1.75...e-02

3000x3000 3.41e-01 3.86...e-02

4000x4000 6.37e-01 7.20...e-02

5000x5000 9.46e-01 1.04...e-01

The slowest is obviously the sliding function since it does 9 samples for ever 1 that block does since I used a 3x3 window...I think...it is the easy stuff that confuses. Enjoy and study the code and its origins should you be interested. I should note that the whole area of sliding implementation has been quite the area of investigation. The implementation that I presented here is one that I chose to use because I liked it. It may not be the fastest or the bestest... but I don't care...if it takes less than a sip of coffee, it is good enough for me.

That's all for now....