Help: How to calculate a fixed specific distance slope?

JulianAlmanzar — Wed, 18 Oct 2023 21:17:08 GMT

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

Re: Help: How to calculate a fixed specific distance slope?

JulianAlmanzar — Wed, 25 Oct 2023 19:49:59 GMT

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)

topic Re: Help: How to calculate a fixed specific distance slope? in ArcGIS Pro Questions

Help: How to calculate a fixed specific distance slope?

Re: Help: How to calculate a fixed specific distance slope?