OREGON STATE UNIVERSITY

Olsen Ex1

Author: Keith Olsen
Email:
Date: 15 March 2011

Here is an example of building a custom moving window analysis with python and arcGIS 10.  The arcGIS 'RasterToNumpyArray' is used to kick a raster into a Numpy array and then the generic_filter in sciPy is used to implement the function.  A custom window shape is used in the example.

# ---------------------------------------------------------------------------
# pixelSearchKernel.py
# Created on: 2011-03-14 15:34:50.00000
#
# Description: This script will run a function in the neighborhood of each pixel
#  of the input raster and output the result as a new raster
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import numpy
from scipy.ndimage import filters

# Local variables:
arcpy.env.workspace = "T:\\Groups\\FPF\\spatial\\IDU\\Deschutes2"
rasterCellSize = 30
rasterInputName = "vegc11_25ac"
rasterOutputName = 'test7'
arcpy.env.cellSize = rasterCellSize
kernelFP = numpy.array([[False,True,False], [True,True,True], [False,True,False]])

# function to perform on kernel
def neighborFunc(kernelArray):
    if (kernelArray[2] != -999):
        return kernelArray[0] + kernelArray[1] + kernelArray[3] + kernelArray[4]
    else:
        return kernelArray[2]

# describe input raster
rasterDesc = arcpy.Describe(rasterInputName)
lowerLeftCorner = arcpy.Point(rasterDesc.Extent.XMin, rasterDesc.Extent.YMin)

# load input raster into numpy array
myArray = arcpy.RasterToNumPyArray(rasterInputName, lowerLeftCorner, rasterDesc.width, rasterDesc.height, -999)

# run neighborFunc over each pixel and save as newRaster
newRaster = filters.generic_filter(myArray, neighborFunc, footprint=kernelFP)

# save numpy array to raster and define projection
rasterToSave = arcpy.NumPyArrayToRaster(newRaster, lowerLeftCorner, rasterCellSize, rasterCellSize, -999)
rasterToSave.save(rasterOutputName)
arcpy.DefineProjection_management(rasterOutputName, rasterDesc.spatialReference)
arcpy.BuildRasterAttributeTable_management(rasterOutputName, "Overwrite")