Creating a Raster Filter

I am looking to create a raster filter which reclassifies single or contiguous pixels less than, for example, 30 pixels (shown in white below) into another class (shown in black). The filter, majority filter, aggregate and other spatial analyst generalization tools seem to be too over reaching in this case. Any Python or Arc solutions and ideas would be appreciated. • Could you be more precise about what you mean by a "cluster"? (Many intuitive and natural definitions are in common use; some differ in important ways.) – whuber Jul 6 '12 at 22:31
• Aaron, what exactly do you mean by "contiguous"? RegionGroup will handle some definitions but will fail for others; morphological operations will succeed for some; and yet other solutions based on focal statistics would be appropriate for other definitions of contiguity. – whuber Jul 9 '12 at 12:22
• In this case, "contiguous" refers to pixels that share a common boundary to the N,S,E,W with a main "group". Pixels adjacent, but diagonal to another pixel would not be considered part of the main group. – Aaron Jul 9 '12 at 12:31
• For that definition, Aaron, @om_henners gives a great solution. – whuber Jul 9 '12 at 12:40

Depending on how complex your raster is, you could run Region Group over the raster (assuming Spatial Analyst here). Then reclassify any region with a Count <= 30 as the black value. Unforunately though Region Group is limited in the number of groups to the size of a table entry an ArcInfo Raster can have.

Alternately you can try the same thing in Python using Numpy and Scipy, and with a lot more control, using the scipy.ndimage.measurements.label function to group the features, then using numpy.bincount to count the values, before reclassifying the orginal raster. Of course this means you have to read your raster into a Numpy array. But assuming you do, code could be as follows:

#first assume our original data is in numpy.array a, with black values as 0, white as 1

#then import our libraries
import numpy as np
from scipy.ndimage.measurements import label

#define how raster cells touch
connection_structure = np.array([[0,1,0],
[1,1,1],
[0,1,0]])

#perform the label operation
labelled_array, num_features = label(a)

#Get the bincount for all the labels
b_count = np.bincount(labelled_array.flat)

#Then use the bincount to set the output values
out_array = np.where(b_count[labelled_array] <= 30, 0, a) #Set small groups to black,
#otherwise maintain original vals

Regarless of how its done though, be prepared as it may take some time to process depending on array complexity.