focal statistics with variable radius

Question

is there an easy way to do focal statistics on a raster, with a variable radius for the neighbourhood to be searched? As in: that search radius would be stored in another raster?

i tried something like

f_max = arcpy.sa.FocalStatistics(Tau,arcpy.sa.NbrCircle(rad_ras,"CELL"), "MAXIMUM","DATA")

it returns

Traceback (most recent call last):
  File "K:\Informatik\tools\arcpy\GKPROZ\test_3.py", line 383, in <module>
    ufereros_gq(r"F:\25sg\25107\TG7\access\gis_logger_tg7.mdb",gq_nummer, False)
  File "K:\Informatik\tools\arcpy\GKPROZ\test_3.py", line 300, in ufereros_gq
    f_max = arcpy.sa.FocalStatistics(Tau,arcpy.sa.NbrCircle(rad_ras,"CELL"), "MAXIMUM","DATA")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4796, in FocalStatistics
    ignore_nodata)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4783, in wrapper
    neighborhood = Utils.compoundParameterToString(neighborhood, ParameterClasses._Neighborhood)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 75, in compoundParameterToString
    return str(parameter)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\ParameterClasses.py", line 202, in __str__
    self.units])
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\CompoundParameter.py", line 39, in _toString
    userProvidedPounds = [value == "#" for value in values]
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 3598, in EqualTo
    in_raster_or_constant2)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 3595, in wrapper
    return _wrapLocalFunctionRaster(u"EqualTo_sa", ["EqualTo", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset # does not exist or is not supported

although the raster(s) exists. I assume this means that the Nbr function can't work with a raster, as the same line works fine for a fixed value for the radius of the circle. I have found a workaround in that i do the focal statistics 10 times with a different radius each time and then use arcpy.sa.Pick to choose from those 10. This is however painfully slow. Is there a better way to do it?

any help would be appreciated.

i use ArcGis 10.0 SP4 on an info license on Windows 7 64

Answer 1

There's no "easy" way to do it that I can see unfortunately. I think you've got the right idea with stacking rasters with the correct focal distance and picking the answer. However Spatial Analyst seems a bit slow...

I wrote a quick script to test the times taken if I did the processing in scipy.ndimage (version 0.7 to be compatible with arcpy) to compare it with Spatial Analyst.

My results were on random arrays of floats were a startling difference between ArcGIS and native Python (even taking into account dumping the raster to a NumPy Array with arcpy):

array size             |    10x10   |   100x100  |   200x200   | 1000x1000
method                 |            |            |             |
ndimage generic filter | 0.540 secs | 4.722 secs | 18.762 secs | FAILED
ndimage max filter     | 0.002 secs | 0.009 secs | 0.034  secs | 0.780 secs
arcpy focal stats and  |            |            |             |
pick                   | 6.744 secs | 3.418 secs | 3.445 secs  | FAILED

My suspicision is for the first time you use a Spatial Analyst tool it loads it into memory for the duration of the script, which is why the time taken drops for the 100x100 and the 200x200 (if I run it outside a loop all three are roughly the same 6.7 seconds).

Run against 1000 points we start to run into limitations with ArcGIS and arcGIS ran out of memory (even with a pre-generated raster on disk).

In the meantime I suspect you're better off using the native maximum filter in scipy.ndimage if you can. That said, to do so you will need to manage memory as depending on the size of your raster and the format you don't want to have to read it all into memory at once!

At any rate an example of sample code for the variable filter:

import numpy as np
import numpy.ma as ma
from scipy import ndimage
import arcpy

a = arcpy.NumPyArrayToRaster(np.random.random((size, size)))
a_mask = arcpy.NumPyArrayToRaster(np.random.randint(0,5, (size, size)))

in_raster = arcpy.RasterToNumPyArray(a)
in_mask = arcpy.RasterToNumPyArray(a_mask)

max_filters = np.dstack([ndimage.filters.maximum_filter(in_raster,
                                                        footprint=get_circle(i, i),
                                                        mode="constant", cval=0)
                         for i in np.unique(in_mask)])
out_arr = np.take(max_filters, in_mask)
arcpy.NumPyArrayToRaster(out_arr).save("somewhere_on_disk.tif")

Answer 2

I think you've got the right idea with stacking rasters with the correct focal distance and picking the answer. However Spatial Analyst seems a bit slow...

I bet you can make it a lot faster by creating a mask from your points using from EucAllocation. Set that mask before you run FocalStatistics so you don't need to do your expensive FocalStatistics run on the entire input raster. You could create one mask for each buffer distance -- but even if you used a single mask generated using your maximum search distance there is probably a lot of time to be saved.