Sliding/Moving window operations in rasters and arrays

07-07-2016 02:16 AM
Labels (1)
MVP Esteemed Contributor
1 0 4,189

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, 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.


import numpy as np
from numpy.lib.stride_tricks import as_strided
from functools import wraps
from textwrap import dedent

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.
    : 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
    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"""

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 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....

About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...