[GIS] Calculating slope cell-by-cell in ArcPy

arcpypython-2.6slope

I am relatively new to the ArcGIS-Python environment and want to calculate elevation raster slope. But it says "ValueError: setting an array element with a sequence". Below code is what I have.

import arcpy, numpy, math
from numpy import *
from arcpy import env

arcpy.env.workspace = ("D:/input")
inRas = arcpy.RasterToNumPyArray('t1')
[rows,cols]=inRas.shape
slope=zeros((rows,cols), float)

desc = arcpy.Describe('t1')
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)

for j in range(0,rows-1):
    for i in range(0,cols-1):
        slope[i,j]= inRas[math.sqrt(((cols+1)-(cols-1)) ** 2 + ((rows+1)-(rows-1)) ** 2)]
        print slope
        myArray=arcpy.NumPyArrayToRaster('slope',pnt,cellSize,cellSize)

        myArray.save("D:/output/t1")

Best Answer

One way to do it is to use the scipy.ndimage generic_filter function. Especially if as you say you're looking at slight variant on the slope function. The following is a quick example of implementing slope using generic_filter on 3x3 footprint, assuming you have the numpy array and the x and y cell sizes (note, you still need to filter for null values and worry about edge effects which are handled in the mode argument of the filter function):

import numpy
from scipy.ndimage import generic_filter

def calc_slope(in_filter, x_cellsize, y_cellsize):
    #slope calculation here - note need to reshape in array to be 3*3
    if -9999 in in_filter:
        return -9999 #Will return -9999 around edge with mode constant and cval -9999
    else:
        [a, b, c, d, e, f, g, h, i] = in_filter
        #From 3x3 box, row 1: a, b, c
        #              row 2: d, e, f
        #              row 3: g, h, i

        dz_dx = ((c + 2*f + i) - (a + 2 * d + g)) / (8 * float(x_cellsize))
        dz_dy = ((g + 2*h + i) - (a + 2 * b + c)) / (8 * float(y_cellsize))
        slope = numpy.sqrt(dz_dx ** 2 + dz_dy**2)

        return numpy.degrees(slope) #we want slope in degrees rather than radians

#slope will return a numpy array the sme size and dimensions as the input raster_data
slope = generic_filter(raster_data, calc_slope, size=3, mode='constant',
                       cval=-9999, extra_arguments=(x_cellsize, y_cellsize))
Related Question