Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have two intersecting line featureclasses. I want to find the angle at each point of intersection using ArcGIS 10 and Python. Can anyone help?

share|improve this question
    
Anyone see anything wrong with this? I keep getting a syntax error.Thanks. math.asin(abs(( !Dx! * !Dya! – !Dy! * !Dxa! ))/( !R! * !Ra! ))/math.pi*180 –  JayB Nov 22 '13 at 22:06
    
If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. –  BradHards Nov 22 '13 at 22:27

2 Answers 2

There is a relatively simple workflow. It overcomes the potential problems that two features may intersect in more than one point. It does not require scripting (but can readily be turned into a script). It can be done primarily from the ArcGIS menu.

The idea is to exploit a layer of intersection points, one point for each distinct pair of intersecting polylines. You need to obtain a small piece of each intersecting polyline at these intersection points. Use the orientations of these pieces to compute their intersection angles.

Here are the steps:

  1. Make sure each of the polyline features has a unique identifier within its attribute table. This will be used later to join some geometric attributes of the polylines to the intersection point table.

  2. Geoprocessing|Intersect obtains the points (make sure to specify you want points for the output).

  3. Geoprocessing|Buffer lets you buffer the points by a tiny amount. Make it really tiny so that the portion of each line within a buffer does not bend.

  4. Geoprocessing|Clip (applied twice) limits the original polyline layers to just the buffers. Because this produces new datasets for its output, subsequent operations will not modify the original data (which is a good thing).

    Figure

    Here is a schematic of what happens: two polyline layers, shown in light blue and light red, have produced dark intersection points. Around those points tiny buffers are shown in yellow. The darker blue and red segments show the results of clipping the original features to these buffers. The rest of the algorithm works with the dark segments. (You cannot see it here, but a tiny red polyline intersects two of the blue lines at a common point, producing what appears to be a buffer around two blue polylines. It's really two buffers around two overlapping points of red-blue intersection. Thus, this diagram displays five buffers in all.)

  5. Use the AddField tool to create four new fields in each of these clipped layers: [X0], [Y0], [X1], and [Y1]. They will hold point coordinates, so make them doubles and give them lots of precision.

  6. Calculate Geometry (invoked by right-clicking on each new field header) enables you to compute the x- and y- coordinates of the start and end points of each clipped polyline: put these into [X0], [Y0], [X1], and [Y1], respectively. This is done for each clipped layer, so 8 calculations are needed.

  7. Use the AddField tool to create a new [Angle] field in the intersection point layer.

  8. Join the clipped tables to the intersection point table based on common object identifiers. (Joins are performed by right-clicking on the layer name and selecting "Joins and Relates".)

    At this point the point intersection table has 9 new fields: two are named [X0], etc., and one is named [Angle]. Alias the [X0], [Y0], [X1], and [Y1] fields which belong to one of the joined tables. Let's call these (say) "X0a", "Y0a", "X1a", and "Y1a".

  9. Use the Field Calculator to compute the angle in the intersection table. Here's a Python code block for the calculation:

    dx = !x1!-!x0!
    dy = !y1!-!y0!
    dxa = !x1a!-!x0a!
    dya = !y1a!-!y0a!
    r = math.sqrt(math.pow(dx,2) + math.pow(dy,2))
    ra = math.sqrt(math.pow(dxa,2) + math.pow(dya,2))
    c = math.asin(abs((dx*dya - dy*dxa))/(r*ra)) / math.pi * 180
    

    The field calculation expression is, of course, merely

    c
    

Despite the length of this code block, the math is simple: (dx,dy) is a direction vector for the first polyline and (dxa,dya) is a direction vector for the second. Their lengths, r and ra (computed via the Pythagorean Theorem), are used to normalize them to unit vectors. (There ought to be no problem with zero lengths, because clipping should produce features of positive length.) The size of their wedge product dx*dya - dy*dxa (after division by r and ra) is the sine of the angle. (Using the wedge product rather than the usual inner product should provide better numerical precision for near-zero angles.) Finally, the angle is converted from radians to degrees. The result will lie between 0 and 90. Note the avoidance of trigonometry until the very end: this approach tends to produce reliable and easily computed results.

Some points may appear multiple times in the intersection layer. If so, they will get multiple angles associated with them.

Buffering and clipping in this solution are relatively expensive (steps 3 and 4): you don't want to do it this way when millions of intersection points are involved. I have recommended it because (a) it simplifies the process of finding two successive points along each polyline within the neighborhood of its intersection point and (b) buffering is so basic it is easy to do in any GIS--no additional licensing is needed above the basic ArcMap level--and usually produces correct results. (Other "geoprocessing" operations might not be so reliable.)

share|improve this answer
    
This could work, but you can't reference field names in the codeblock, so you'd have to wrap the code in a function and call it using the field names as arguments. –  mvexel Apr 21 '12 at 5:52
    
@mv Thank you for that observation. One could also use VBS instead of Python--VBS will parse field names in the code block. –  whuber Apr 22 '12 at 20:38
1  
It actually worked like a charm when using a function wrapper. I found that in ArcGIS 10 and when using Python, you don't need to alias the variables, you can prepend the join table name in the field reference, like !table1.x0!. –  mvexel Apr 24 '12 at 16:39

I believe you need to create python script.

You can do it using geoprocessing tools and arcpy.

Here is the main tools and ideas that can be useful for you:

  1. Make intersection of your two polyline (lets call them PLINE_FC1, PLINE_FC2) featureclasses (you need point features as a result - POINT_FC) using tool Intersect. You will have IDs from PLINE_FC1, PLINE_FC2 in points POINT_FC.
  2. Split PLINE_FC1 by POINT_FC using tool Split Line At Point. By result you will have splitted polylines - main advantage of it that you can just take first/last vertex of such line compare it to the next/previous vertex (coordinate difference) and calculate angle. So, you will have angle of your line in the point of intersect. There is one problem here - you have to run this tool manually several times to realize how output is written. I mean if it takes polyline, split it, write two result polylines to the output and then proceed to the next polyline and repeat. Or may be this part (result of splitting) are written to different memory feature classes, then appended to the output. This is the main problem - to realize how output is written to be able to filter only first part of each polyline after splitting. Another possible solution is to loop throught all result splitted polylines with SearchCursor and take only first encountered (by ID of source polylines PLINE_FC1).
  3. In order to get angle you will need to access result polyline's verteces using arcpy. Write resulting angles to points (POINT_FC).
  4. Repeat steps 2-3 for PLINE_FC2.
  5. Substract angle attributes (in POINT_FC) and get result.

May be it will be very difficult to code step 2 (also some tools require ArcInfo license). Then you can also try to analyse verteces of every polyline (grouping them by ID after intersection).

Here is the way to do it:

  1. Take first intersection point POINT_FC. Get its coordinates (point_x, point_y)
  2. By its ID take respective source polyline from PLINE_FC1.
  3. Take first (vert0_x, vert0_y) and second (vert1_x, vert1_y) verteces of it.
  4. For the first vertex calculate tangent of line between this vertex and intersection point: tan0 = (point_y - vert0_y) / (point_x - vert0_x)
  5. Calculate same thing for second vertex: tan1 = (vert1_y - point_y) / (vert1_x - point_x)
  6. If tan1 is equal to tan2, then you have found two verteces of your line which have intersection point in between and you can calculate angle of intersection for this line. Otherwise you have to proceed to the next pair of verteces (second, third) and so on.
  7. Repeat steps 1-6 for every intersection point.
  8. Repeat steps 1-7 for second polyline featureclass PLINE_FC2.
  9. Substract angle attributes from PLINE_FC1 and PLINE_FC2 and get result.
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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