Filter... Convolve... Rolling... Sliding

DanPatterson_Retired · ‎11-26-2017

Filter

Started this with rolling stats.

Rolling Statistics for grid data... an alternative

Once you have a filter... kernel whatever you want to call it, all you need is a moving window of some form and shape. I won't go into those details since my previous work covers this. I will just report a few simple ones, and show you some of the possible kernels that you can create.

# ---- my input array ----

Input array.......
-shape (1, 10, 10), ndim 3
.... 1 1 1 1 1 2 2 2 2 2  
.... 1 1 1 1 1 2 2 2 2 2  
.... 1 1 1 1 1 2 2 2 2 2  
.... 1 1 1 1 1 2 2 2 2 2  
.... 1 1 1 1 1 2 2 2 2 2  
.... 3 3 3 3 3 4 4 4 4 4  
.... 3 3 3 3 3 4 4 4 4 4  
.... 3 3 3 3 3 4 4 4 4 4  
.... 3 3 3 3 3 4 4 4 4 4  
.... 3 3 3 3 3 4 4 4 4 4‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

From the above, you can see that I kept it simple... 4 blocks of values, the perfect study area.

Let's see if we can pick up some lines that divide the raster/array

Sobel horizonal filter...
-shape (1, 3, 3), ndim 3

  . -1 -2 -1  
  .  0  0  0  
  .  1  2  1  


Sobel horizontal result.......
-shape (1, 10, 10), ndim 3, masked array, fill value 999999

.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1  8  8  8  8  8  8  8  8 -1  
.... -1  8  8  8  8  8  8  8  8 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1  0  0  0  0  0  0  0  0 -1  
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

See those 8's! nailed the horizontal divide.

vertical filter...

  . -1  0  1  
  . -2  0  2  
  . -1  0  1  

Sobel vertical result.......
-shape (1, 10, 10), ndim 3, masked array, fill value 999999
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1  0  0  0  4  4  0  0  0 -1  
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So easy... vertical nailed as well.

One more...

Array...
-shape (1, 3, 3), ndim 3
  . -1 -1  0  
  . -1  0  1  
  .  0  1  1  

Emboss vertical result.......

-shape (1, 10, 10), ndim 3, masked array, fill value 999999
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1  4  4  4  6  6  4  4  4 -1  
.... -1  4  4  4  6  6  4  4  4 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1  0  0  0  2  2  0  0  0 -1  
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

How about some filters that I have found.

    n = np.nan
    all_f = [1, 1, 1, 1, 1, 1, 1, 1, 1]
    no_cnt = [1, 1, 1, 1, n, 1, 1, 1, 1]
    cross_f = [1, 0, 1, 0, 0, 0, 1, 0, 1]
    plus_f = [0, 1, 0, 1, 0, 1, 0, 1, 0]
    grad_e = [1, 0, -1, 2, 0, -2, 1, 0, -1]
    grad_n = [-1, -2, -1, 0, 0, 0, 1, 2, 1]
    grad_ne = [0, -1, -2, 1, 0, -1, 2, 1, 0]
    grad_nw = [2, -1, 0, -1, 0, 1, 0, 1, 2]
    grad_s = [1, 2, 1, 0, 0, 0, -1, -2, -1]
    grad_w = [-1, 0, 1, -2, 0, 2, -1, 0, 1]
    lap_33 = [0, -1, 0, -1, 4, -1, 0, -1, 0]
    line_h = [-1, -1, -1, 2, 2, 2, -1, -1, -1]
    line_ld = [2, -1, -1, -1, 2, -1, -1, -1, 2]
    line_rd = [-1, -1, 2, -1, 2, -1, 2, -1, -1]
    line_v = [-1, 0, -1, -1, 2, -1, -1, 2, -1]
    high = [-0.7, -1.0, -0.7, -1.0, 6.8, -1.0, -0.7, -1.0, -0.7]
    sob_hor = [-1, -2, -1, 0, 0, 0, 1, 2, 1]   # sobel y
    sob_vert = [-1, 0, 1, -2, 0, 2, -1, 0, 1]  # sobel x
    emboss = [-1, -1, 0, -1, 0, 1, 0, 1, 1]
    sharp1 = [0., -0.25, 0., -0.25, 2.0, -0.25, 0., -0.25, 0.]
    sharp2 = [-0.25, -0.25, -0.25, -0.25, 3.0, -0.25, -0.25, -0.25, -0.25]
    sharp3 = [-1, -1, -1, -1, 9, -1, -1, -1, -1]
    lowpass = [1, 2, 1, 2, 4, 2, 1, 2, 1]

    # ---- assemble the dictionary ----

    d = {1: all_f, 2: no_cnt, 3: cross_f, 4: plus_f, 5: grad_e,
         6: grad_n, 7: grad_ne, 8: grad_nw, 9: grad_s, 10: grad_w,
         11: lap_33, 12: line_h, 13: line_ld, 14: line_rd, 15: line_v,
         16: high, 17: sob_hor, 18: sob_vert, 19: emboss, 20: sharp1,
         21: sharp2, 23: sharp3, 24: lowpass}

    filter_ = np.array(d[mode]).reshape(3, 3)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

You call the filter using the mode number (ie d[23] is sharp3.

import numpy as np
from numpy.lib.stride_tricks import as_strided

def stride(a, win=(3, 3), stepby=(1, 1)):
    """Provide a 2D sliding/moving view of an array.
    :  There is no edge correction for outputs. Use _pad first.
    """
    err = """Array shape, window and/or step size error.
    Use win=(3,3) with stepby=(1,1) for 2D array
    or win=(1,3,3) with stepby=(1,1,1) for 3D
    ----    a.ndim != len(win) != len(stepby) ----
    """
    assert (a.ndim == len(win)) and (len(win) == len(stepby)), err
    shape = np.array(a.shape)  # array shape (r, c) or (d, r, c)
    win_shp = np.array(win)    # window      (3, 3) or (1, 3, 3)
    ss = np.array(stepby)      # step by     (1, 1) or (1, 1, 1)
    newshape = tuple(((shape - win_shp) // ss) + 1)
    newshape += tuple(win_shp)
    newstrides = tuple(np.array(a.strides) * ss) + a.strides
    a_s = as_strided(a, shape=newshape, strides=newstrides).squeeze()
    return a_s


def a_filter(a, mode=1, ignore_nodata=True, nodata=None):
    """put the filter stuff below this line"""
    #
    # ---- stride the input array ----
    a_strided = stride(a)
    if ignore_nodata:
        c = (np.sum(a_strided * filter_, axis=(2, 3)))
    else:
        c = (np.nansum(a_strided * filter_, axis=(2, 3)))
    pad_ = nodata
    if nodata is None:
        if c.dtype.name in ('int', 'int32', 'int64'):
            pad_ = min([0, -1, a.min()-1])
        else:
            pad_ = min([0.0, -1.0, a.min()-1])
    c = np.lib.pad(c, (1, 1), "constant", constant_values=(pad_, pad_))
    m = np.where(c == pad_, 1, 0)
    c = np.ma.array(c, mask=m, fill_value=None)
    return filter_, c 

# ----------------------------------------------------------------------
# __main__ .... code section
if __name__ == "__main__":
    """Optionally...
    : - print the script source name.
    : - run the _demo
    """
#    you can put your stuff here... use RasterToNumPyArray to get an array‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So that is pretty well some other stuff.

I should point out that the same process can be used for statistical calculations, terrain derivatives (slope, aspect, hillshade, curvature etc ) and I have reported on these before.

Now I should point out that none of the above requires the Spatial Analyst extension and the RasterToNumPyArray and NumPyArrayToRaster are provided with arcpy.da 's module. Lots of possibilities, lots of capability is built-in.

More later.

SteveKopp · ‎11-29-2017

Thanks Dan, nice examples. Didn't see it in your earlier post so just wanted to point out that if anyone is reading this and doesn't want to write code for edge detection and other filter types, the easiest approach is to use the Convolution raster function which has 20 built in filter types and a UI for creating custom filters in the dialog.

And for those who prefer writing python, there are some nice Jupyter notebook examples for raster, including this one I quite like which shows how to use scikit image to count center pivot irrigation fields from an image. For python devs looking to go beyond ArcGIS image analysis capabilities scikit image is a good package to explore.

DanPatterson_Retired · ‎11-29-2017

Thanks Steve, good sources to be sure. I left out the 'pad' option for the blogs since the method of padding would depend on the filter being used. The simplest is of course to just put a nodata collar around the raster.

The convolution function's builtins are definitely a source of some of the kernels being shown. Before they became 'free', they required the spatial analyst extension, so some of this is historic.

As for scikit image... yes, but not for ArcMap types, with Pro it is much easier to work with the Sci stack.

Too bad we can't have Jupyter notebooks in our blogs... play and show... Jive isn't up to speed on that yet. For educational purposes, I would love to see notebooks part of the experience (posted on this around UC.

SteveKopp · ‎11-29-2017

Yes embedded Jupyter would be great, and would save people like you a lot of formatting hassle.