Generalization tools for rasters...

10-08-2018 03:14 AM
Labels (1)
MVP Esteemed Contributor
0 0 797

Scipy ndimage morphology … I am betting that is the first thought that popped into your head.  A new set of tools for the raster arsenal.

But wait!  Their terminology is different than that used by ArcGIS Pro or even ArcMap.

What tools are there?  Lots of filtering tools, but some of the interesting ones are below

from scipy import ndimage as nd
dir(nd)[... snip ...
'affine_transform', ... , 'center_of_mass', 'convolve', 'convolve1d', 'correlate',
'correlate1d', ... 'distance_transform_edt', ... 'filters', ... 'generic_filter',
... 'geometric_transform', ... 'histogram', 'imread', 'interpolation', ... 'label',
'labeled_comprehension', ... 'measurements', ... 'morphology', ...  'rotate', 'shift',

I previously covered distance_transform_edt which performs the equivalent of Euclidean Distance and Euclidean Allocation in this post


Let's begin with a raster constructed from 3x3 windows with a repeating pattern and repeats of the pattern.

You can construct your own rasters using numpy, then apply a NumPyArrayToRaster to it.  I have covered this before, but here is the code to produce the raster

def _make_arr_():
    """Make an array with repeating patterns

    a = np.arange(1, 10).reshape(3, 3)
    aa = np.repeat(np.repeat(a, 3, axis=1), 3, axis=0)
    aa = np.tile(aa, 3)
    aa = np.vstack((aa, aa))
    return aa‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

---- The raster (image) ----

---- (1) Expand ----

Expand the zone 1 values by 2 distance units

---- (2) Shrink ----

Shrink the zone 1 values by 1 distance unit.

---- (3) Regiongroup ----

Produce unique regions from the zones.  observe that the zone 1 values in the original are given unique values first. Since there were 6 zone 1s, they are numbered from 1 to 6.  Zone 2 gets numbered from 7 to 12.  and that pattern of renumbering repeats.

---- (4) Aggregate ----

Aggregation of the input array using the mode produces a spatial aggregation retaining the original values.  Our original 3x3 kernels are now reduced to a 1x1 size which has 3 times the width and height (9x area).

---- (5) Some functions to perform the above ----

In these examples buff_dist is referring to 'cells'.

For the aggregation, there are a number of options as shown in the header. 

Generally integer data should be restricted to the mode, min, max, range and central (value) since the median and mean upscale to floating point values.  This of course can be accommodated by using the python statistics package median_low and median_high functions:

statistics — Mathematical statistics functions — Python 3.7.1rc1 documentation 

So think of a function that you want.  Filtering is a snap since you can 'stride' an array using any kernel you want using plain numpy or using the builtin functions from scipy.

Enough for now...

def expand_(a, val=1, mask_vals=0, buff_dist=1):
    """Expand/buffer a raster by a distance

    if isinstance(val, (list, tuple)):
        m = np.isin(a, val, invert=True).astype('int')
        m = np.where(a==val, 0, 1)
    dist, idx = nd.distance_transform_edt(m, return_distances=True,
    alloc = a[tuple(idx)]
    a0 = np.where(dist<=buff_dist, alloc, a)  #0)
    return a0

def shrink_(a, val=1, mask_vals=0, buff_dist=1):
    """Expand/buffer a raster by a distance

    if isinstance(val, (list, tuple)):
        m = np.isin(a, val, invert=False).astype('int')
        m = np.where(a==val, 1, 0)
    dist, idx = nd.distance_transform_edt(m, return_distances=True,
    alloc = a[tuple(idx)]
    m = np.logical_and(dist>0, dist<=buff_dist)
    a0 = np.where(m, alloc, a)  #0)
    return a0

def regions_(a, cross=True):
    """currently testing regiongroup
    np.unique will return values in ascending order


    if (a.ndim != 2) or (a.dtype.kind != 'i'):
        msg = "\nA 2D array of integers is required, you provided\n{}"
        return a
    if cross:
        struct = np.array([[0,1,0], [1,1,1], [0,1,0]])
        struct = np.array([[1,1,1], [1,1,1], [1,1,1]])
    u = np.unique(a)
    out = np.zeros_like(a, dtype=a.dtype)
    details = []
    is_first = True
    for i in u:
        z = np.where(a==i, 1, 0)
        s, n = nd.label(z, structure=struct)
        details.append([i, n])
        m = np.logical_and(out==0, s!=0)
        if is_first:
            out = np.where(m, s, out)
            is_first = False
            n_ = n
            out = np.where(m, s+n_, out)
            n_ += n
    details = np.array(details)
    details = np.c_[(details, np.cumsum(details[:,1]))]
    return out, details

def aggreg_(a, win=(3,3), agg_type='mode'):
    """Aggregate an array using a specified window size and an aggregation

    a : array
        2D array to perform the aggregation
    win : tuple/list
        the shape of the window to construct the blocks
    agg_type : string aggregation type
        max, mean, median, min, mode, range, sum, central


    blocks = block(a, win=win)
    out = []
    for bl in blocks:
        b_arr = []
        for b in bl:
            if agg_type == 'mode':
                uni = np.unique(b).tolist()
                cnts = [len(np.nonzero(b==u)[0]) for u in uni]
                idx = cnts.index(max(cnts))
            elif agg_type == 'max':
            elif agg_type == 'mean':
            elif agg_type == 'median':
            elif agg_type == 'min':
            elif agg_type == 'range':
                b_arr.append((np.nanmax(b) - np.nanmin(b)))
            elif agg_type == 'sum':
            elif agg_type == 'central':
                b_arr.append(b.shape[0]//2, b.shape[1]//2)
                tweet("\naggregation type not found {}".format(agg_type))
                b_arr.append(b.shape[0]//2, b.shape[1]//2)
    out = np.array(out)
    return out
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...