Calculating Distance and Direction using ArcPy

/in ,

While working with geospatial information, it is often advantageous to find out how close one particular piece of data is to other pieces of data. This leads to a greater understanding of the area of study. The knowledge of how things relate to one another spatially is articulated in Waldo Tobler’s First Law of Geography. It states that “everything is related to everything else, but near things are more related than distant things.”

This became especially important in one of our current contracts, as we needed to discover which feature in a large infrastructure dataset was closest to a single area of interest. This task can easily be done in ArcMap manually by loading in both datasets, zooming to the area of interest, and using the ruler tool to measure the distance between the area of interest and the nearest point in the infrastructure dataset.

It can also be done by running the near analysis tool in ArcToolbox, using the area of interest and the infrastructure dataset as inputs. However, I wanted to automate the task in order to avoid any human error that may be present in the process.

I decided to create a stand-alone script using python that would both calculate the distance and direction from the area of interest to the nearest piece of infrastructure data. The script required relatively minimal coding in order to achieve the desired output. First, I defined the area of interest and the infrastructure dataset input and output variables:

input1 = r"D:\Path_To...\Area_of_Interest.shp"
input2 = r"D:\Path_To...\Infrastructure_dataset.shp"
output1 = r"D:\Path_To...\Area_of_Interest_Projected.shp"
output2 = r"D:\Path_To...\Infrastructure_dataset_Projected.shp"

I then projected both of the datasets to the same projection, in order to calculate the distance correctly. I chose NAD 1983 UTM Zone 18N in this example, reflected by the WKID 26918, which uses meters as its unit of measurement. Then I ran the near analysis:

arcpy.Project_management(input1, output1, arcpy.SpatialReference(26918))
arcpy.Project_management(input2, output2, arcpy.SpatialReference(26918))
arcpy.Near_analysis(output1, output2, "", "LOCATION", "ANGLE")

Next, I obtained the distance and direction values returned by the near analysis, which added two fields to the projected area of interest dataset. Then I converted the distance from meters to miles:

field1 = "NEAR_DIST"
field2 = "NEAR_ANGLE"
cursor = arcpy.SearchCursor(output1)
for row in cursor:
    distance = (row.getValue(field1))
    angle = (row.getValue(field2))
distance = (distance / 1609.34)

After that, I rounded the distance and angle values to two decimal places and converted the angle from a -180 – 180 degree value to a 0-360 degree value. I chose to convert the angle in order to more easily calculate the direction the nearest infrastructure point is from the area of interest data point.

distance = str(round(distance,2))
angle = round(angle,2)
if angle < 0:
    angle = angle +360
    print angle
else:
    angle = angle

Here I ran an if statement to update the direction variable, based upon the angle value obtained earlier, and printed the distance and direction for the user to see:

direction = ""
if angle >= 0 and angle < 22.5:
    direction = "east"
elif angle >= 22.5 and angle < 67.5:
    direction = "northeast"
elif angle >= 67.5 and angle < 112.5:
    direction = "north"
elif angle >= 112.5 and angle < 157.5:
    direction = "northwest"
elif angle >= 157.5 and angle < 202.5:
    direction = "west"
elif angle >= 202.5 and angle < 247.5:
    direction = "southwest"
elif angle >= 247.5 and angle < 292.5:
    direction = "south"
elif angle >= 292.5 and angle < 337.5:
    direction = "southeast"
elif angle >= 337.5 and angle < 360:
    direction = "east"
print "The direction value returned is "+direction+ "."   
dist_ang = "The nearest piece of infrastructure is approximately "+ distance+ " miles\n"+direction+" of the area of interest."
print dist_ang

I later included and expanded this functionality in another script I had developed and associated it with an ArcGIS tool, which output the resultant calculations into a text element in the map layout.

The full source code for the script discussed in this post can be found on GitHub.

This post was written by:
Christian
Associate

For more information on this post, or our geospatial information services, please e-mail us at contact@zekiah.com.