Slope between max and min values of a defined area within a DEM

Say I have a DEM raster. I then also I have a set of polygons in which I would like to know the slope between the max and min DEM values of any given polygon. If my thinking is correct, I could calculate the max slope as follows...

from within any given polygon:

(maxElev - minElev) / distance(maxElevPoint, minElevPoint)

I think the math is correct, but what tools within 3D Analyst could help me to automate this. I can Python and know ArcPy, but I'm not sure where to start. I found the Slope tool that deals with rasters within the 3D Analyst extension, I could convert to raster but from what I can tell this gives the slope back for each cell compared to it's adjacent cells.

what is the best workflow for doing this type of analysis?

• When trying to figure out a new workflow I always think the best place to start is with some test code. Work from what you know to what you don't and ask about where you get stuck.
– PolyGeo
Jun 21, 2017 at 20:32
• Your math is wrong. What makes you think it give you max slope. Slope raster is a way to go Jun 21, 2017 at 23:31
• @FelixIP Can you please explain how the math is wrong? Jun 23, 2017 at 15:05
• I am not interested in finding the slope between adjacent cells. I.E. I don't want a raster output in which each cell is given the max slope value of its adjacent cells. I want the maximum slope of the portion of the raster contained by each polygon. Jun 23, 2017 at 15:09
• Slope 1 / 1 is bigger than 50/100, i.e. largest difference in elevation does not always produce greatest slope Jun 23, 2017 at 19:18

I truly struggle to understand your definition of slope and believe slope derived drom DEM is the one you need. Alternatively you might apply TIN and take a slope of steepest triangle.

However If you insist on your definition of slope, considering 2 extreme points to derive the highest slope is wrong. Let's assume for simplicity that you have only 4 points sitting on a straight line:

Isn't it obvious, that highest slope is not between points 1,4 (your extremes!), but between points 2 and 3? Thus you'll need to iterate between all combinations of points' pairs. Good luck with this, knowing density of LiDAR points...

SOLUTION BASED ON UPDATED QUESTION, assumes points have PGON_ID assigned:

``````arcpy.Sort_management("LIDAR_POINTS", "../high_points.shp", "Z DESCENDING")
arcpy.DeleteIdentical_management("high_points", "PGON_ID")
arcpy.Sort_management("LIDAR_POINTS", "../low_points.shp", "Z ASCENDING")
arcpy.DeleteIdentical_management("low_points", "PGON_ID")
arcpy.Merge_management("high_points;low_points", "../pairs.shp")
arcpy.PointsToLine_management("pairs","../lines.shp","PGON_ID")
``````

RESULT shows original points labelled by Z:

Note: it is very fast but won't handle duplicated extreme points, it will pick first ones.

• I definitely agree that points 2 and 3 have a higher slope than 1 and 4. I believe you are correct in that using the word "maximum" is not appropriate. I suppose what I'm really after is the slope between the maximum elevation and the minimum elevation for a give area. When/if there are multiple identical max/min values within the area, I want to choose the pair with the furthest distance. I apologize if I wasn't clear. Thank you for the explanation above. Jul 10, 2017 at 11:56
• I think I will change the word "maximum" in the OP. Jul 10, 2017 at 12:00
• Now it makes sense Jul 10, 2017 at 20:22

After a while of messing around with this I came up with a workflow that works. It is not the fastest and I'm sure could be improved. I don't want to accept this answer yet in case someone has anything to add. The general workflow is the following:

For each input polygon:

• extract the portion of the raster that falls beneath the polygon
• convert the extracted raster to points
• get the min and max value of all points
• find the min/max pair with the greatest distance between them
• calculate slope, angle of incline, max difference, etc.
• insert a new polygon that contains all of the new stat data

Here is the fully commented script:

``````import arcpy, os, math

#Input tool parameters
dem = arcpy.GetParameterAsText(1)
outPath = arcpy.GetParameterAsText(2)
outName = arcpy.GetParameterAsText(3)

#create empty fc for the results to be stored in
arcpy.CreateFeatureclass_management(outPath, outName, "POLYGON", "", "", "", arcpy.Describe(inPads).spatialReference)

srDEM = arcpy.Describe(dem).spatialReference

if sr.linearUnitName != srDEM.linearUnitName:
arcpy.AddError("Linear Units of Input Spatial References Do Not Match! If the linear units don't fit, I must aquit!!!  Goodbye...")
else:
arcpy.AddWarning("Spatial Reference Check Passed...  You may proceed!")

#add necessary fields to the new empty FC

#SQL query constructor
def buildWhereClauseFromList(table, field, valueList):
"""Takes a list of values and constructs a SQL WHERE
clause to select those values within a given field and table."""

# Determine field type
fieldType = arcpy.ListFields(table, field)[0].type

# Add single-quotes for string field values
if str(fieldType) == 'String':
valueList = ["'%s'" % value for value in valueList]

# Format WHERE clause in the form of an IN statement
whereClause = "%s IN(%s)" % (fieldDelimited, ', '.join(map(str, valueList)))
return whereClause

#list for data that needs to be eliminated after the tool runs.
deleteList = []

#turn the input polys into an "in memory layer"

#Set a count int
id = 0

#Basic arcpy.da search cursor implementation to cycle through each feature in the original poly fc
with arcpy.da.SearchCursor('pads_lyr', ["OBJECTID", "SHAPE@"]) as cursor:
for row in cursor:

#construction of a where clause to actually select each row for further geo-processing

#clear current selection - prob not necessary, but why not...

#make a row selection based on the current OID used to construct the where clause (SQL)

#copy the selected poly feature to a temporary fc that contains only obe feature.

#clip the portion of the dem that falls under the geometry of the current row
arcpy.Clip_management(dem, "", os.path.join(outPath, 'tempRast' + str(id)), os.path.join(outPath, 'tempPad' + str(id)), "", "ClippingGeometry")

#clear current selection - again, prob not necessary...

#convert the raster cells to points with values
arcpy.RasterToPoint_conversion(os.path.join(outPath, 'tempRast' + str(id)), os.path.join(outPath, 'tempPoints' + str(id)))

#Get XY

#Get min and max values from point features
min = arcpy.SearchCursor(os.path.join(outPath, 'tempPoints' + str(id)), "", "", "", "grid_code" + " A").next().getValue("grid_code")
max = arcpy.SearchCursor(os.path.join(outPath, 'tempPoints' + str(id)), "", "", "", "grid_code" + " D").next().getValue("grid_code")

arcpy.MakeFeatureLayer_management(os.path.join(outPath, 'tempPoints' + str(id)), 'points_lyr')

whereMax = buildWhereClauseFromList('points_lyr', "grid_code", [max])
whereMin = buildWhereClauseFromList('points_lyr', "grid_code", [min])

maxCoords = []
minCoords = []

minPoint = "NULL"
maxPoint = "NULL"

#start tie breaking distance, this gets recycled later as the result becomes the greatest excepted distance
dist = 0.0

#Select all of the points with the maxValue and append their coordinates to the coordinate list
arcpy.SelectLayerByAttribute_management('points_lyr', "NEW_SELECTION", whereMax)
with arcpy.da.SearchCursor('points_lyr', ["POINT_X", "POINT_Y"]) as tCurs:
for r in tCurs:
maxCoords.append([r[0], r[1]])

#Select all of the points with the minValue and append their coordinates to the coordinate list
arcpy.SelectLayerByAttribute_management('points_lyr', "NEW_SELECTION", whereMin)
with arcpy.da.SearchCursor('points_lyr', ["POINT_X", "POINT_Y"]) as tCurs:
for r in tCurs:
minCoords.append([r[0], r[1]])

#In case of multiple points with the same min or max values, use the largest distanced coordinates
for maxCoord in maxCoords:
maxGeom = arcpy.PointGeometry(arcpy.Point(maxCoord[0], maxCoord[1]), sr)
for minCoord in minCoords:
minGeom = arcpy.PointGeometry(arcpy.Point(minCoord[0], minCoord[1]), sr)

if minGeom.distanceTo(maxGeom) > dist:
dist = minGeom.distanceTo(maxGeom)
minPoint = minGeom
maxPoint = maxGeom

arcpy.AddMessage("Excepted Coordinates: " + "MaxValuePoint(" + str(maxPoint.firstPoint.X) + ", " + str(maxPoint.firstPoint.Y) + ")  MinValuePoint(" + str(minPoint.firstPoint.X) + ", " + str(minPoint.firstPoint.Y) + ")")

#get the difference of the min and max..  Min and Max never change even if a further distanced pair is selected from tie breaking
diff = float(max - min)

#get slope value for max and min points
slope = float(diff/dist)

#get angle of incliniation ----> Using simple trig inverser tangent, we know that tan(angle) = rise/run.
angle = math.degrees(math.atan((diff/dist)))

#construct a "data access" inserter for the new result FC and all the necessary fields
inserter = arcpy.da.InsertCursor(os.path.join(outPath, outName), ["SHAPE@", "MAX_ELEV", "MIN_ELEV" ,"MAX_CHANGE", "MAX_X_COORD", "MAX_Y_COORD", "MIN_X_COORD", "MIN_Y_COORD", "DIST", "MAX_SLOPE", "INCLN_ANGLE", "ORIG_OID"])

#use the new inserter to insert a new row ( note:  the shape@ token gets passed directly from the original search cursor...  is this a no, no, i.e. nested cursors? - john davies)
inserter.insertRow([row[1], max, min, diff, maxPoint.firstPoint.X, maxPoint.firstPoint.Y, minPoint.firstPoint.X, minPoint.firstPoint.Y, dist, slope, angle, row[0]])

#delete the inserter from RAM
del inserter

arcpy.AddMessage("                                         max/min elevation difference: " + str(diff))
arcpy.AddMessage("                                                     max/min distance: " + str(dist))

#append the delete list with all temp data from this iteration
deleteList.append(os.path.join(outPath, 'tempPoints' + str(id)))
deleteList.append(os.path.join(outPath, 'tempRast' + str(id)))
#delete "in mem" layer
arcpy.Delete_management('points_lyr')

#up that counter
id = id + 1

arcpy.AddMessage("All Features have been Calculated...  Gotta clean up a bit... this might take a second, or two, or three...  maybe a lot more...")

#join the original poly FCs data back to the new result FC (in case there are any fields that should be preserved)

#cleanup everything in the delete list
for item in deleteList:
arcpy.Delete_management(item)

#cleanup some uneccesary fields
arcpy.DeleteField_management(os.path.join(outPath, outName), "ORIG_OID")

#cleanup more uneccesary fields
try:
arcpy.DeleteField_management(os.path.join(outPath, outName), "SHAPE_Length_1")
except:
#WHy cant delete these fields???
arcpy.AddWarning("An expected cleanup field 'SHAPE_Length_1' wasn't found...  Everything is probably okay, but check your results.")

try:
arcpy.DeleteField_management(os.path.join(outPath, outName), "SHAPE_Area_1")
except:
arcpy.AddWarning("An expected cleanup field 'SHAPE_Length_1' wasn't found...  Everything is probably okay, but check your results.")

``````

This works really well. Again, the downfall is how long it takes to run if there are many polygons.

• I'd do a `raster to point` and then a `spatial join` between the points and the polygons for the entire raster and polygon. Don't do it for each feature. Jul 10, 2017 at 15:28