Help: How to calculate a fixed specific distance slope?

[JulianAlmanzar]

I have a DEM with a horizontal resolution of 1x1m. I need to calculate the slope, but instead of the usual slope calculation with neighboring cells, I require the slope to be determined as follows:

• Each cell will calculate the slope with all other cells at exactly a 20m distance.
• Only the maximum slope from these calculations is considered.

I have been unable to specify a fixed distance for slope calculation, so I attempted to write a code to perform this process (I converted the raster to points beforehand):

File1="points"
File2="points_1"
Field="slope"

with arcpy.da.SearchCursor(File2, ["SHAPE@", "OBJECTID","grid_code"]) as cursor:

    for row in cursor:

        objectid = row[1]
        sql = "OBJECTID = {}".format(objectid)
        arcpy.management.SelectLayerByAttribute(File2, "NEW_SELECTION", sql)
        val=row[2]

        arcpy.management.SelectLayerByLocation(File1, 'WITHIN_A_DISTANCE', File2,'19.5 Meters', 'NEW_SELECTION')
        arcpy.management.SelectLayerByAttribute(File1, 'SWITCH_SELECTION')
        arcpy.management.SelectLayerByLocation(File1, 'WITHIN_A_DISTANCE', File2,'20.5 Meters', 'SUBSET_SELECTION')

        with arcpy.da.SearchCursor(File1, "grid_code") as cursor2:
            values=[]
            for row2 in cursor2:
                value=row2[0]
                values.append(value)
            max_value=max(values)
            min_value=min(values)

        slope=max(abs(max_value - val),abs(min_value - val))/5

        arcpy.management.CalculateField(File2, Field, slope)

Things to note:

• "points" and "points_1" are identical copies. Having two copies was the only way I found to be able to select the points at 20m (actually between 19.5m and 20.5m).
• I am aware that I still need to use trigonometric functions to calculate the slope in degrees.

The code works as intended; however, the number of points from almost any region is colossal. I attempted to extract just a small area to test the code, but that 'small' area resulted in 523568 points, which means the code would take weeks to process them all. I need to be able to do this for much larger areas, so it is currently unviable.

I need help to achieve slope calculations without taking up a significant amount of time. I appreciate any assistance. Thanks.

Julian

[JulianAlmanzar]

After many tests, I finally found a way to do it faster. I think it is still not as efficient as it could be, but now it processes the half million points in around 20 min, instead of several months. The solution I got was using numpy and arrays instead of the selection tools from arcgis, here is how it looks:

import arcpy

import numpy as np

import math

radius = 20

# Transform the raster into an NumPy array

raster = "LiDAR"

raster_array = arcpy.RasterToNumPyArray(raster)

num_rows, num_cols = raster_array.shape

# Extracts the spatial coordinates for the previously created array

spatial_reference = arcpy.Describe(raster).spatialReference

cell_size_x = arcpy.Describe(raster).meanCellWidth

cell_size_y = arcpy.Describe(raster).meanCellHeight

x_min = arcpy.Describe(raster).extent.XMin

y_max = arcpy.Describe(raster).extent.YMax

y_min = arcpy.Describe(raster).extent.YMin

# Create two new arrays for the X and Y coordinates each

x_coords = np.arange(x_min, x_min + num_cols * cell_size_x, cell_size_x)

y_coords = np.arange(y_max, y_max - num_rows * cell_size_y, -cell_size_y)

CoordX, CoordY = np.meshgrid(x_coords, y_coords)

# calculate the slope

slope_array = np.empty_like(raster_array) # creates an empty array for the slope

for row in range(num_rows) :

    for col in range(num_cols) :

        value = raster_array[row,col] # Extracts the elevetation value

        Value_X = CoordX[row,col]   # Gets the X Y coordinates

        Value_Y = CoordY[row,col]

        row_start= max(0, row - radius)

        row_end= min(num_rows, row + radius + 1)

        col_start = max(0, col - radius)

        col_end = min(num_cols, col + radius + 1)

        # Loop to find the cells at an specific radius distance, and obtain the max and min values of elevation from those cells

        minvalue = 1000000

        maxvalue = 0

        row_start= max(0, row - radius) # limiting search area to increase speed

        row_end= min(num_rows, row + radius + 1)

        col_start = max(0, col - radius)

        col_end = min(num_cols, col + radius + 1)

        for row2 in range(row_start,row_end) :

            for col2 in range(col_start,col_end) :

                distance = math.sqrt((CoordX[row2,col2] - Value_X)**2 + (CoordY[row2,col2] - Value_Y)**2)

                if (radius-0.5) <= distance <= (radius+0.5) :

                    val=raster_array[row2,col2]

                    if val > maxvalue:

                        maxvalue = val

                    if val < minvalue:

                        minvalue = val

        #slope calculation to degrees

        slope=max(abs(maxvalue - value),abs(minvalue - value))/radius

        slopes_degrees = np.degrees(np.arctan(slope))

        slope_array[row,col]=slopes_degrees # updates the slope array

slope_raster = arcpy.NumPyArrayToRaster(slope_array, arcpy.Point(x_min, y_min), cell_size_x, cell_size_y)