Can the Particle Track (Spatial Analyst) tool be terminated when a specific change in the path direction angle occurs?

4271
7
02-19-2015 07:41 AM
MatthewPinson
New Contributor

I developed a Python script to read through a CSV file of coordinate points and use those coordinate points in the Particle Track (Spatial Analyst) tool. Now that I have created the particle track polylines associated with the list of coordinate points, I need to truncate the polylines if the particle track path changes direction by more than 90 degrees. I looked through all of the available toolsets, but I cannot identify one that will trim the polylines in this manner.

 

I thought of a couple of possible solutions, but I do not know how to implement them (or if it is even possible). One method would be to somehow customize the Particle Track tool to end its process if the track direction changes by more than 90 degrees. Alternately, I considered somehow adjusting the input raster to create an artificial zone of zero magnitude that would stop the tracking process, but I have not figured out how I can specify which raster cells to alter.

 

Any suggestions?

0 Kudos
7 Replies
DarrenWiens2
MVP Honored Contributor

Not exactly sure I've captured what you want, but this script will cycle through all the lines in a polyline feature class, and output lines up until it reaches an acute angle.

>>> import math
... lines = []
... with arcpy.da.SearchCursor("YOUR_POLYLINE_FC_HERE","SHAPE@") as cursor:
...    for row in cursor:
...        for part in row[0]:
...            linePnts = []
...            writePnt = 1
...            pt_count = 1
...            for pnt in part:
...                pnt = arcpy.PointGeometry(pnt)
...                if pt_count > 1:
...                    if pt_count >2:
...                        distAB = oneBack.distanceTo(twoBack)
...                        distBC = pnt.distanceTo(oneBack)
...                        distAC = pnt.distanceTo(twoBack)
...                        angB = math.degrees(math.acos((((distAB*distAB)+(distBC*distBC))-(distAC*distAC))/(2*distAB*distBC)))
...                        if angB < 90:
...                            writePnt = 0
...                    twoBack = oneBack
...                oneBack = pnt
...                pt_count += 1
...                if writePnt == 1:
...                    linePnts.append(arcpy.Point(pnt.centroid.X,pnt.centroid.Y))
...            line = arcpy.Polyline(arcpy.Array(linePnts))
...            lines.append(line)
... arcpy.CopyFeatures_management(lines,r'in_memory\lines')

1.PNG

0 Kudos
MatthewPinson
New Contributor

Thank you for the suggestion, but I don't know what this code means or how I can use it with my existing data. I honestly know very little about Python, or programming in general, so some comments or description that explicitly define and instruct how to use the code might be helpful. I don't know if I could just use the code as is by substituting in the name of the feature class file that contains the polylines, or if I would need to customize the code in some way to suit my needs. I see the indication of the substitution of the feature class file, but what is meant by "SHAPE@"? What is the output of this code and where is it saved? Does it overwrite the existing feature class? If it is necessary to customize the code for my particular situation, I will definitely need some additional guidance. Thank you.

0 Kudos
DarrenWiens2
MVP Honored Contributor

In ArcMap, copy/paste the code into the Python window (under Geoprocessing menu). Change the name of the layer in the SearchCursor (keep within quotes). It will output to an in_memory layer called lines and add it to your map. You can export that layer to one on disk, or change the script to save to a new location (e.g. substitute r'C:\someFolder\someShapefile.shp' for r'in_memory\lines'). "SHAPE@" is a geometry token - don't change it. You can think of it as "the geometry".

0 Kudos
MatthewPinson
New Contributor

No luck. The new lines created by the code appear exactly the same as the existing lines. See the screen grab. The new lines (ultra blue) overlay the original lines (malachite green). I made the line style wider for the original lines so they can be seen beneath the new lines, but otherwise they appear to trace the same path. Thanks for trying though.

line_trim_test.jpg

0 Kudos
DarrenWiens2
MVP Honored Contributor

Hmmm not sure why it doesn't work for you. Is your line feature class projected? If not, that may affect the calculated distances between points, which would affect the angle calculations.

0 Kudos
MatthewPinson
New Contributor

The line feature class is a merge of the individual polylines created by the Particle Track tool. The code finally worked to trim the lines after I ran the Unsplit Lines tool; however, it stripped the feature class attributes, and, in some instances, left some disconnected line segments floating in open space. I can easily disregard the floating line segments with a select by location that only includes those lines that intersect the boundary area. It would be helpful to be able to trim the lines without losing the feature class attributes, though.

0 Kudos
MatthewPinson
New Contributor

Any ideas why the code might encounter a math domain error? At this point, I've created new lines with updated input data and I am trying to test the code for trimming the lines on this new feature class. As before, the code does an exact trace of the lines that were directly output by the Particle Track tool (without catching the angle change of significance). My earlier assumption was that the executed code did not find the acute angles due to the segmented format of the lines drawn by the Particle Track tool (essentially each line segment was treated independently, thus no acute angles existed in the process of the code). I believed that I had accommodated for this by running the Unsplit Line tool before executing the tracing/trimming code. That method seemed to work with my previous run of the data; however, using the tracing/trimming code on the unsplit version of the new lines generates a math domain error (apparently in some part of the "angB = math.degrees(math.acos..." calculation).

0 Kudos