# Calculating slope cell-by-cell in ArcPy?

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")
``````
-
Look at the line where you calculate slope and notice that your formula does not depend on the loop indices `i` and `j`! (The formula simplifies to `inRas[math.sqrt(8)]`.) Indeed, just what does it mean to index a raster (`inRas`) by a floating point value? – whuber Dec 30 '11 at 15:33
If you wanted to loop across the entire array your could try ndenumerate instead of the nested for loops... – om_henners Jan 2 '12 at 0:27

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))
``````
-
As a suggestion, I would "float" the x_cellsize and y_cellsize or change the "8" in the dz_dx and dz_dy equations to "8.0" to avoid integer division should the cell sizes be passed as integers in error or by omission. – Dan Patterson Dec 31 '11 at 14:48
Cheers @DanPatterson - edited to float the x and y cell sizes. – om_henners Jan 2 '12 at 0:22
thanks om_henners it will surely help. – Ibe Jan 3 '12 at 16:10
``````    slope[i,j]= inRas[math.sqrt(((cols+1)-(cols-1)) ** 2 + ((rows+1)-(rows-1)) ** 2)]
``````

I believe that you probably wanted something like:

``````    slope[i,j]= math.sqrt(((inRas[j+1])-(inRas[j-1])) ** 2 + ((inRas[i+1])-(inRas[i-1])) ** 2)]
``````

Note that I read that real quick, so it could easily be off (i/j swapped, or similar), but as far as I can tell you are trying to basically do sqrt(a^2 + b^2) where a and b are the differences between horizontally and vertically adjacent pixel values centered on the target pixel. If that's right, then something closer to the above is probably what you intended. Hope that helps.

-
That's closer, but it still won't fly, for at least two reasons. (1) When i or j are 0, `inRas[j-1]` or `inRas[i-1]` will be out-of-bounds references. (2) The formula takes no account of cellsize. This is a nice example of why one should prefer using the built-in capabilities of the GIS rather than rolling one's own: ArcGIS has a capable, fast, accurate slope calculation that can be invoked in a single line of code. – whuber Dec 30 '11 at 16:43
Absolutely correct. Technically, you have to inset the result by 1/2 the kernel size in both directions for any of these masking types of operations (and adjust appropriately for even/odd kernel sizes). I was simply dealing with the basic array referencing issue. That said, I agree with you. Don't re-invent the wheel. If the function exists and can be used for your purposes, use it. – Brian Dec 30 '11 at 19:26
thanks. The only reason to reinvent wheel is that i have to write a slight variant of slope like which is different than hard coded in arcgis. Therefore above snippet was a starter. And as pointed the code is giving an index error due to out of range processing. Can you please hint at how to inset? – Ibe Dec 30 '11 at 23:30
Limiting the range of the loops from 1..cols-1 and 1..rows-1 will take care of the "inset". You should first initialize the slopes to null values (not zeros) so that the border of cells where slopes are not calculated is effectively clipped out. – whuber Dec 31 '11 at 14:26
thanks whuber for pointing to this issue. – Ibe Jan 3 '12 at 16:11