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.
You must be a registered user to add a comment. If you've already registered, sign in. Otherwise, register and sign in.