Calculation of directional rotation angles from 3 points

Question

What I'm trying to do is find the angle of rotation between lines formed by three consecutive points. These are sequential points, so the direction of rotation matter. My input is a sequence of pair of coordinates.

My desired output is the rotation angle of each point, where the point would be acting as the vertex of the angle. This angle would be between 1 and 360, where a negative number indicates a rotation to the left, and a positive number a rotation to the right.

I've been struggling with this for weeks, but I'm finally closer to the solution than ever before. I wrote the following script and compared it to the output of the "pathmetrics" function of the program Geospatial Modelling Tool(GME).

coords=coords=[[283907,971700],[284185,971634],[284287,971507],[284275,971608],[283919,971761],[284311,971648],[284277,971637],[284280,971651],[284174,971649],[283909,971701],[283941,971700],[284294,971518],[284288,971517],[284315,971539],[284250,971505]
print "A"+"\t"+"B"+"\t"+"C"+"\t"+"orientation"+"\t"+"angle"+"\t"+"bearing AB"+"\t"+"bearing BC"
for a in coords:
  A=a
  indexA=coords.index(a)
  B=coords[indexA+1]
  C=coords[indexA+2]
  ##Find the bearings of AB and BC
  AB=[B[0]-A[0],B[1]-A[1]]          #find the extreme of vector AB
  BearAB=math.atan2(AB[0],AB[1])    #use arctan2 to find the angle
  ABBearDeg=math.degrees(BearAB)    #in degrees
  if ABBearDeg<0:                   #if negative, add 360 in order to obtain the angle in a clockwise direction
   ABBearDeg=360+ABBearDeg          #Bearing AB
  BC=[C[0]-B[0],C[1]-B[1]]          #Do the same for points BC
  BearBC=math.atan2(BC[0],BC[1])
  BCBearDeg=math.degrees(BearBC)
  if BCBearDeg<0:
   BCBearDeg=360+BCBearDeg          #Bearing BC
 ##Find the angle between the lines
  alfa=BCBearDeg-ABBearDeg          #Obtain the difference between the bearing angles
  if abs(alfa)>180:                 #If greater than 180
   if alfa<0:                        #and negative
    angle=(360+alfa)                   #Direction of rotation is right and angle is obtained by adding 360
    print format(A)+"\t"+format(B)+"\t"+format(C)+"\t"+"right"+"\t"+format(angle)+"\t"+format(round(ABBearDeg,2))+"\t"+format(round(BCBearDeg,2))
   else:                             #If positive
    angle=alfa-360                      #Direction of rotation is left and angle is obtained by substracting 360
    print format(A)+"\t"+format(B)+"\t"+format(C)+"\t"+"left"+"\t"+format(angle)+"\t"+format(ABBearDeg)+"\t"+format(round(BCBearDeg,2))
  else:                            #If the difference was less than 180, then the rotation angle is equal to it
   angle=alfa
   if angle<0:                     #If negative, left rotation
       print format(A)+"\t"+format(B)+"\t"+format(C)+"\t"+"left"+"\t"+format(angle)+"\t"+format(ABBearDeg)+"\t"+format(round(BCBearDeg,2))
   else:                            #If positive, right rotation
    print format(A)+"\t"+format(B)+"\t"+format(C)+"\t"+"right"+"\t"+format(angle)+"\t"+format(ABBearDeg)+"\t"+format(round(BCBearDeg,2))

While many of the results coincide, others don't.

                            My results                                                                      GME results     
        A                   B                   C       orientation     angle       bearing AB  bearing BC  GMEangle    GMEbearing AB   GMEbearing BC
[283907, 971700]    [284185, 971634]    [284287, 971507]    right       37.8750     103.3553    141.2300    37.8750     103.3553        141.2303
[284185, 971634]    [284287, 971507]    [284275, 971608]    left        -148.0060   141.2303    353.2200    -148.0060   141.2303        353.2243
[284287, 971507]    [284275, 971608]    [283919, 971761]    left        -59.9675    353.2243    293.2600    -68.9673    353.2243        284.2570
[284275, 971608]    [283919, 971761]    [284311, 971648]    right       172.8236    293.2600    106.0800    96.6181     284.2570        106.0804
[283919, 971761]    [284311, 971648]    [284277, 971637]    right       145.9916    106.0804    252.0700    145.9916    106.0804        252.0721
[284311, 971648]    [284277, 971637]    [284280, 971651]    right       120.0227    252.0700    12.0900     120.0227    252.0721        12.0948
[284277, 971637]    [284280, 971651]    [284174, 971649]    left        -103.1757   12.0948     268.9200    -103.1757   12.0948         268.9191
[284280, 971651]    [284174, 971649]    [283909, 971701]    right       12.1828     268.9191    281.1000    -131.4097   268.9191        137.5094
[284174, 971649]    [283909, 971701]    [283941, 971700]    right       170.6880    281.1000    91.7900     85.2053     137.5094        9.4623
[283909, 971701]    [283941, 971700]    [284294, 971518]    right       25.4848     91.7899     117.2700    146.4722    9.4623          85.0429
[283941, 971700]    [284294, 971518]    [284288, 971517]    right       143.2629    117.2748    260.5400    123.0283    85.0429         260.5377
[284294, 971518]    [284288, 971517]    [284315, 971539]    right       150.2887    260.5400    50.8300     150.2887    260.5377        50.8263
[284288, 971517]    [284315, 971539]    [284250, 971505]    left        -168.4394   50.8263     242.3900    -147.5925   50.8263         263.2338

(sorry for the messy table; I couldn't upload an image as I wanted)

I've been able to pin-point at with point the error occurs, but I can't figure out why it happens because it depends strictly on pre-set formulas I have no control of. So, the difference is that (SOMETIMES) my calculations of the bearing of the vectors differs from the one calculated by GME. The weird part is that it happens only sometimes, and I have no idea what triggers it.

Any ideas of what might be going on?

If you know any other way of calculating angles between lines that incorporate the direction of the movement, do let me know. Anything that works will be okay with me.

Thanks!!!!

this may help http://stackoverflow.com/questions/13226038/calculating-angle-between-two-lines-in-python — Joran Beasley, Mar 15 '13 at 21:04

Roland Smith · Answer 1 · 2013-03-15T22:00:15.697

You have the points A, B and C. I'm assuming from your description that B is the rotation point. That is, the vector BA transforms into BC.

Create the vectors BA (A-B) and BC (C-B).
Calculate the angle of the vectors (make them positive) and their difference
The difference in angles between the vectors is the rotation angle.

Using numpy makes this easy. Like this:

In [1]: import numpy as np

In [2]: A = np.array([283907, 971700])

In [3]: B = np.array([284185, 971634])

In [4]: C = np.array([284287, 971507])

In [5]: BA = A - B

In [6]: BC = C - B

In [7]: BA
Out[7]: array([-278,   66])

In [8]: BC
Out[8]: array([ 102, -127])

In [9]: s = np.arctan2(*BA)

In [10]: if s < 0:
   ....:     s += 2*np.pi
   ....:     

In [11]: s
Out[11]: 4.9454836529138948

In [12]: e = np.arctan2(*BC)

In [13]: if e < 0:
    e += 2*np.pi
   ....:     

In [14]: e
Out[14]: 2.4649341681747883

In [15]: delta = e - s

In [16]: np.degrees(delta)
Out[16]: -142.12501634890182

In [17]: delta
Out[17]: -2.4805494847391065

A positive angle is counter-clockwise.

Calculation of directional rotation angles from 3 points

1 Answers1

Linked