nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Calculating angles using ArcGIS Desktop 10.2

BrianVann

I need to find angles of streets where they intersect with roadways covering the entire state. Does anyone have a method and/or tools in mind?

300px-Angle_obtuse_acute_straight.svg.png

タグ付けされた投稿をもっと見る

Accepted answers

WesMiller

Attached is a model that should accomplish what you are trying to do.If you'll add your streets layer and leave the defaults it should produce a point file at each intersection that has an angle between 5 and 40. I made these parameters so they could be adjusted, the other distance field is how far back from the intersection do you want to check the angle, 50 was used you may want smaller. The model will need an Advanced license to run.

findAngle.tbx.zip

All comments

SumvegMalik

Hi Wes

I was trying to use the toolbox but I'm not able to use it after some operations. Snapshot is attached below:

it is not going forward after find less than angle. also could you explain what is Angle lower and Angle upper in the parameters

Thanks in advance

Sumveg

WesMiller

You could run Dissolve—Help | ArcGIS for Desktop then Spatial Join—Help | ArcGIS for Desktop either can be used at all license levels. When you do your spatial join you'll need to work with Mapping input fields to output fields—Help | ArcGIS for Desktop to join your street names

JeremyMullins4

Just as an amendment to Darren's awesome script, I noticed that my output was including intersections at endpoints of all line segments, and not just between streets. For example, in my feature class, Main St was made up of ten different line segments (for E911 purposes).

To counter act this, I simply added in the following Dissolve tool before finding the intersection (changes in bold😞

out_fc = r'C:\Users\jmullins\Desktop\GEO\GEOImportTool.gdb\Scratch\points'

int_pt = r'C:\Users\jmullins\Desktop\GEO\GEOImportTool.gdb\Scratch\int_pt'

diss_fc = r'C:\Users\jmullins\Desktop\GEO\GEOImportTool.gdb\Scratch\diss_fc'

arcpy.Dissolve_management(fc,diss_fc,["YourStreetNameField"],"","SINGLE_PART")

arcpy.Intersect_analysis(diss_fc,int_pt,output_type='POINT')

This will dissolve all streets by their street name, which should decrease the number of redundant or unnecessary points in your output. Many may not use Darren's script for this in particular, especially if they aren't interested in street intersections, but I just wanted to add it in there for those who may need it.

Thanks again, everyone.

JeremyMullins4

Wes,

I mentioned in my comment that I tried that and it fixed the first error. However, the second error came up afterwords. I tried changing the paths from 'in_memory' to a file geodatabase and it gave me the same error. I looked through the script provided by ArcGIS for Intersect in the location provided in the error, and I didn't see any syntax issues (obviously, why would there be if it worked for Darren and others...)

I see where in line 289 it gives me an error, but I'm not sure why. Here is what the script says:

from arcpy.geoprocessing._base import gp, gp_fixargs

from arcpy.arcobjects.arcobjectconversion import convertArcObjectToPythonObject

try:

retval = convertArcObjectToPythonObject(gp.SpatialJoin_analysis(*gp_fixargs((target_features, join_features, out_feature_class, join_operation, join_type, field_mapping, match_option, search_radius, distance_field_name), True)))

return retval

except Exception, e:

raise e

I'm trying to figure out why it raised e. I am only using a single polyline feature class to perform this analysis, but I didn't see that as an issue with the Intersect tool. I made sure to replace my lines to the path to the file geodatabase (although I have a feeling it may need to be a shapefile to access the .dbf file???)

WesMiller

Put the below at the top of the script

import arcpy

JeremyMullins4

Darren,

I'm getting the following issue running your python script:

Traceback (most recent call last):

File "C:\Users\jmullins\Desktop\Test.py", line 3, in <module>

sr = arcpy.Describe(fc).spatialReference

NameError: name 'arcpy' is not defined

I threw an import arcpy underneath import os, and now I am getting this error:

Traceback (most recent call last):

File "C:\Users\jmullins\Desktop\Test.py", line 8, in <module>

arcpy.Intersect_analysis(fc,int_pt,output_type='POINT')

File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\analysis.py", line 289, in Intersect

raise e

RuntimeError: Object: Error in executing tool

I apologize for being inept at python scripting, and I appreciate any help you (or anyone else) can provide. Thanks!

BrianVann

Darren,

I finally got the script to work by exporting the fc into a shapefile.

And now receiving this error message:

Runtime error Traceback (most recent call last): File "<string>", line 49, in <module> File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py",

line 3052, in AddField raise e ExecuteError: ERROR 000732: Input Table: Dataset in_memory\points does not exist or is not supported

Brian

BrianVann

Darren,

I also thought of completing this task similar to what your image shows. It seems it should be somewhat simple; buffer points at intersections by 100 ft, convert the buffers to a line, and simply clip the arc of each resulting line where the roads intersect the circular line. Any circular line that has been clipped from a roadway line that has a arc length between the parameters of 1ft to 122.17ft (at a radius of 100ft) would be an intersection where a roadway has an intersecting roadway of less than 70 degrees. The problem I had with this method is that the process returned many arcs lengths from the buffered lines that were actually not due to the clipping function but rather due to the fact that they were not completely closed lines and thus left short arcs...

Your model looks very promising and seems a good method to use; however, I keep getting error messages when put into the python window. I haven’t much experience with python.

example: Runtime error Traceback (most recent call last): File "<string>", line 9, in <module> File "c:\program files (x86) \arcgis\desktop10.2\arcpy\arcpy\management.py", line 4021, in Dissolve raise e ExecuteError: ERROR 000369: Invalid input field(s)

Hope this makes sense...

And, thanks for the help!!

DarrenWiens2

The script makes points on each buffer in isolation, so there is no influence from any other buffer. I agree that the buffer size should be as low as reasonable. You would run into problems if the buffer overlapped other intersections or parallel streets.

Good:

Bad:

WesMiller

Darren NICE!

Does your second intersect create points where the buffers overlap see above from Brian? I told Brian to reduce the buffer size and that should work, but like you I've thought a lot about this one too. In retrospect i think i would add a select by location and just gather the points that intersect the streets and calculate from there. After looking at your code i feel like i got off easy using the near tool to grab the angles of the line and doing light math to calculate the angles.

DarrenWiens2

This post has been on my mind for a while, and here's what I came up with. I suspect it's quite similar to Wes' previously posted solution, but works for all license levels:

>>> import os
... fc = "myLines"
... sr = arcpy.Describe(fc).spatialReference
... radius = 50
... out_fc = r'in_memory\points'
... int_pt = r'in_memory\int_pt'
... arcpy.Intersect_analysis(fc,int_pt,output_type='POINT')
... diss_int_pt = r'in_memory\diss_int_pt'
... arcpy.Dissolve_management(int_pt,diss_int_pt,'#',[["FID","MIN"]],"SINGLE_PART")
... buff = r'in_memory\buff'
... arcpy.Buffer_analysis(diss_int_pt,buff,str(radius) + ' METERS')
... buff_line_int = r'in_memory\buff_line_int'
... arcpy.Intersect_analysis([buff,fc],buff_line_int,output_type='POINT')
... sing_buff_line_int = r'in_memory\sing_buff_line_int'
... arcpy.MultipartToSinglepart_management(buff_line_int,sing_buff_line_int)
... new_points = {}
... with arcpy.da.SearchCursor(diss_int_pt,['OID@','SHAPE@'],spatial_reference=sr) as cursor1:
...    for row1 in cursor1:
...        cent_pt = row1[1].centroid
...        angs = []
...        with arcpy.da.SearchCursor(sing_buff_line_int,'SHAPE@','\"FID_buff\" = ' + str(row1[0]),spatial_reference=sr) as cursor2:
...            for row2 in cursor2:
...                buff_pt = row2[0].centroid
...                dx = cent_pt.X - buff_pt.X
...                dy = cent_pt.Y - buff_pt.Y
...                if dx < 0 and dy <= 0:
...                    ang = math.degrees(math.atan(abs(dy/dx)))
...                if dx <= 0 and dy > 0:
...                    ang = math.degrees(math.atan(abs(dx/dy))) + 270
...                if dx > 0 and dy >= 0:
...                    ang = math.degrees(math.atan(abs(dy/dx))) + 180
...                if dx >= 0 and dy < 0:
...                    ang = math.degrees(math.atan(abs(dx/dy))) + 90
...                angs.append(ang)
...        angs.sort()
...        for i in range(1,len(angs)):
...            mid_ang = ((angs + angs[i-1])/2)
...            ang_diff = angs - angs[i-1]
...            new_x = cent_pt.X + (radius * math.cos(math.radians(mid_ang)))
...            new_y = cent_pt.Y + (radius * math.sin(math.radians(mid_ang)))
...            new_point = arcpy.PointGeometry(arcpy.Point(new_x, new_y),sr)
...            new_points[str(row1[0]) + '_' + str(i)] = [ang_diff,new_point]
...        mid_ang = (((360-angs[-1]) + angs[0])/2) - (360-angs[-1])
...        ang_diff = (360-angs[-1]) + angs[0]
...        new_x = cent_pt.X + (radius * math.cos(math.radians(mid_ang)))
...        new_y = cent_pt.Y + (radius * math.sin(math.radians(mid_ang)))
...        new_point = arcpy.PointGeometry(arcpy.Point(new_x, new_y),sr)
...        new_points[str(row1[0]) + '_0'] = [ang_diff,new_point]
... arcpy.CreateFeatureclass_management(os.path.dirname(out_fc),os.path.basename(out_fc),'POINT',spatial_reference=sr)
... arcpy.AddField_management(out_fc,'FID_buff',"LONG")
... arcpy.AddField_management(out_fc,'ANGLE',"DOUBLE")
... iCursor = arcpy.da.InsertCursor(out_fc,['SHAPE@','FID_buff','ANGLE'])
... for k,v in new_points.iteritems():
...    row = [v[1],k.split('_')[0],v[0]]
...    iCursor.insertRow(row)

WesMiller

Would you be willing to share a portion of your data? 22 is far less than i expected for a state wide search. I'll check the code against your data. I tested this code against county data and returned over 200.

BrianVann

Wes,

Thanks!

I ran the model and it failed with this error:

Assembling Features...
Succeeded at Mon Jan 25 11:06:26 2016 (Elapsed Time: 2 minutes 4 seconds)
Executing (Generate Near Table): GenerateNearTable in_memory\roadsIntersect in_memory\Roads_Intersect_Buffer_Inter in_memory\nearTable "51 Feet" NO_LOCATION ANGLE ALL 0 PLANAR
Start Time: Mon Jan 25 11:06:27 2016
Determining data processing extents...
Building a neighborhood index from the Near Features...
Generating Near Table...
Found 346 feature(s) within 51.000000
Succeeded at Mon Jan 25 11:06:27 2016 (Elapsed Time: 0.56 seconds)
Executing (findLessThanAngle): findLessThanAngle 50 5 in_memory\roadsIntersect in_memory\nearTable
Start Time: Mon Jan 25 11:06:27 2016
Running script findLessThanAngle...
Failed script findLessThanAngle...

Traceback (most recent call last):
File "Y:\GISData\Tools\FindAngletoolbox\findAngle.tbx#findLessThanAngle.py", line 37, in <module>
NameError: name 'each' is not defined

Failed to execute (findLessThanAngle).
Failed at Mon Jan 25 11:06:28 2016 (Elapsed Time: 0.23 seconds)

I deleted the "each" out of the python code, ran it, but it only gave me 22 intersecting roadways with angles of < than 50 degrees out of 2299 total...I wasn't sure if it was my python code editing that gave me the short list.

I know very little about python. It is in my toolbox of things to do, no pun intended.

Other than that, it does look promising!

Thanks again for your help!

WesMiller

findAngle.tbx.zip

DanPatterson_Retired

You might want want to move this to Geoprocessing GeoNet Community Structure this space is for help on GeoNet and how to use it.

BrianVann

Thanks for the help Wes!

It looks as though this tool returns angle values in relation to North, South, East, or West. I am looking for a method that returns an angle where 2 lines intersect. In addition, I have a least a few hundreds of thousands of intersecting lines to measure those angles... We are trying to find all the intersecting roadways in the state that are at an angle of less than 40 degrees, where they intersect, and run crash rates and analyses. Do you still think this tool will do the job of measuring those angles?

Thanks again for the prompt response!

WesMiller

if you have an Advanced License you could use the Near—Help | ArcGIS for Desktop

会話に参加する

サインインしたメンバーは投稿、更新のフォローなどができます。初めてですか？無料アカウントを登録してください。

クイックリンク

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.