0

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 '17 at 20:32
  • Your math is wrong. What makes you think it give you max slope. Slope raster is a way to go – FelixIP Jun 21 '17 at 23:31
  • @FelixIP Can you please explain how the math is wrong? – GeoJohn Jun 23 '17 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. – GeoJohn Jun 23 '17 at 15:09
  • Slope 1 / 1 is bigger than 50/100, i.e. largest difference in elevation does not always produce greatest slope – FelixIP Jun 23 '17 at 19:18
2

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:

enter image description here

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:

enter image description here

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. – GeoJohn Jul 10 '17 at 11:56
  • I think I will change the word "maximum" in the OP. – GeoJohn Jul 10 '17 at 12:00
  • Now it makes sense – FelixIP Jul 10 '17 at 20:22
0

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
inPads = arcpy.GetParameterAsText(0)
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)

sr = 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
arcpy.AddField_management(os.path.join(outPath, outName), "MAX_X_COORD","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MAX_Y_COORD","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MIN_X_COORD","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MIN_Y_COORD","DOUBLE")

arcpy.AddField_management(os.path.join(outPath, outName), "MAX_ELEV","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MIN_ELEV","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MAX_CHANGE","DOUBLE")

arcpy.AddField_management(os.path.join(outPath, outName), "DIST","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "INCLN_ANGLE","DOUBLE")
arcpy.AddField_management(os.path.join(outPath, outName), "MAX_SLOPE","DOUBLE")

arcpy.AddField_management(os.path.join(outPath, outName), "ORIG_OID","LONG")

arcpy.AddMessage("Max Elevation Change Calculation Starting For: " + str(inPads))
arcpy.AddMessage(" ")

#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."""

    # Add DBMS-specific field delimiters
    fieldDelimited = arcpy.AddFieldDelimiters(arcpy.Describe(table).path, field)

    # 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"
arcpy.MakeFeatureLayer_management(inPads, 'pads_lyr')

arcpy.AddMessage(" ")
arcpy.AddMessage(" ")

#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
        where = buildWhereClauseFromList('pads_lyr', "OBJECTID", [row[0]])

        #clear current selection - prob not necessary, but why not...
        arcpy.SelectLayerByAttribute_management('pads_lyr', "CLEAR_SELECTION")

        #make a row selection based on the current OID used to construct the where clause (SQL)
        arcpy.SelectLayerByAttribute_management('pads_lyr', "NEW_SELECTION", where)

        #copy the selected poly feature to a temporary fc that contains only obe feature.
        arcpy.CopyFeatures_management('pads_lyr', os.path.join(outPath, 'tempPad' + str(id)))

        #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...
        arcpy.SelectLayerByAttribute_management('pads_lyr', "CLEAR_SELECTION")

    #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
    arcpy.AddXY_management(os.path.join(outPath, 'tempPoints' + str(id)))

    #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(" ")
        arcpy.AddMessage("Feature Added ID: " + str(id))
        arcpy.AddMessage("                                         max/min elevation difference: " + str(diff))
        arcpy.AddMessage("                                                     max/min distance: " + str(dist))
        arcpy.AddMessage("                                                                slope: " + str(slope))
        arcpy.AddMessage("                                                              incline: " + str(angle))


        arcpy.AddMessage(" ")

    #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)))
        deleteList.append(os.path.join(outPath, 'tempPad' + str(id)))
    #delete "in mem" layer
    arcpy.Delete_management('points_lyr')

    #up that counter
    id = id + 1

arcpy.AddMessage(" ")
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...")
arcpy.AddMessage(" ")
arcpy.AddMessage("...Cleaning Up")
arcpy.AddMessage(" ")

#join the original poly FCs data back to the new result FC (in case there are any fields that should be preserved)
arcpy.JoinField_management(os.path.join(outPath, outName), "ORIG_OID", inPads, "OBJECTID")

#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:
    arcpy.AddMessage(" ")
    #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.AddMessage(" ")
    arcpy.AddWarning("An expected cleanup field 'SHAPE_Length_1' wasn't found...  Everything is probably okay, but check your results.")

arcpy.AddMessage("All Done!")
arcpy.AddMessage(" ")

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. – Emil Brundage Jul 10 '17 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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