Solved! Go to Solution.
def getCardinal(angle): lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"], [157.5, "NW"], [202.5, "W"], [247.5, "SW"], [292.5, "S"], [337.5, "SE"], [360, "E"]]
def getCardinal(angle): lstCardinal = [[11.25,"E"], [33.75,"ENE"], [56.25, "NE"], [78.75, "NNE"], [101.25, "N"], [123.75, "NNW"], etc... [348.75, "ESE"], [360, "E"]]
import arcpy from math import atan2, pi def getCardinal(angle): lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"], [157.5, "NW"], [202.5, "W"], [247.5, "SW"], [292.5, "S"], [337.5, "SE"], [360, "E"]] for item in lstCardinal: value = item[0] if angle < value: cardinal = item[1] break return cardinal def calcGeomCardinality(polyline): pnt1 = polyline.firstPoint pnt2 = polyline.lastPoint angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi if angle_deg < 0: angle_deg = 360 + angle_deg return getCardinal(angle_deg)
calcGeomCardinality( !SHAPE!)
import arcpy from math import atan2, pi def getAngle(polyline): pnt1 = polyline.firstPoint pnt2 = polyline.lastPoint angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi if angle_deg < 0: angle_deg = 360 + angle_deg return angle_deg
getAngle( !SHAPE!)
Hi Mark,
Although this type of analysis is normally done using the raster format, it can be done using a featureclass (shapefile). Assuming that you are interested in an angle between the start point and the end point of each feature you could do this:
- open the attribute table of your shapefile
- add a new field that will hold the cardinal flow direction, make sure it is a text field (let's call this field "CardFD")
- right click on your new field "CardFD" and select "Field Calculator..." to open the Field Calculator
- Choose Python as parser
- Switch the "Show Codeblock" on
- Paste the code below in the "Pre-Logic Script Code":
import arcpy from math import atan2, pi def getCardinal(angle): lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"], [157.5, "NW"], [202.5, "W"], [247.5, "SW"], [292.5, "S"], [337.5, "SE"], [360, "E"]] for item in lstCardinal: value = item[0] if angle < value: cardinal = item[1] break return cardinal def calcGeomCardinality(polyline): pnt1 = polyline.firstPoint pnt2 = polyline.lastPoint angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi if angle_deg < 0: angle_deg = 360 + angle_deg return getCardinal(angle_deg)
- Paste the line below in as formula (just below "CardFD =")
calcGeomCardinality( !SHAPE!)This will take each polyline feature and extract the from and to node from it. Next it will calculate the angle and parse it to cardinal flow direction.
- Click "OK"
You can also write the angle itself to a new field:
- add a new field (double or float) that can hold the angle between the start and end point . Let's call this field "FlowAngle"
- right click on your new field "FlowAngle" and select "Field Calculator..." to open the Field Calculator
- it will probably show the code from cardinal flow direction. Replace the code in the "Pre-Logic Script Code" by
import arcpy from math import atan2, pi def getAngle(polyline): pnt1 = polyline.firstPoint pnt2 = polyline.lastPoint angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi if angle_deg < 0: angle_deg = 360 + angle_deg return angle_deg
- Paste the line below in as formula (just below "FlowAngle=")
getAngle( !SHAPE!)
- Click "OK"
If you are going to perform some analysis to determine the predominant angle, please be aware that the resulting angle does not have to be representative for the area.
Kind regards,
Xander
import arcpy def getSLDist(geom): # Get the end points of the polyline fp = geom.firstPoint tp = geom.lastPoint # Create an Array object and add points to it array = arcpy.Array() array.add(fp) array.add(tp) # Create a polyline from array, as it only has 2 points in it # it must be a straight line, then return its length line = arcpy.Polyline(array) dist = line.length return dist
Mark,
You can compute the straight line distance of any polyline with the following python field calculate. First create a field of type double then run a calculate with the python expression getSLDist(!shape!) and the pre-logic script code is:import arcpy def getSLDist(geom): # Get the end points of the polyline fp = geom.firstPoint tp = geom.lastPoint # Create an Array object and add points to it array = arcpy.Array() array.add(fp) array.add(tp) # Create a polyline from array, as it only has 2 points in it # it must be a straight line, then return its length line = arcpy.Polyline(array) dist = line.length return dist
Duncan