5

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")
Improve this question
2
  • 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, 2011 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, 2012 at 0:27

2 Answers 2

Reset to default
6

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))
Improve this answer
2
  • 2
    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.
    – user681
    Dec 31, 2011 at 14:48
  • Cheers @DanPatterson - edited to float the x and y cell sizes.
    – om_henners
    Jan 2, 2012 at 0:22
1 
    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.

Improve this answer
4
  • 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, 2011 at 16:43
  • 2
    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, 2011 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, 2011 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, 2011 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.