### Thread: Rotation angle using 3 points

1. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Mar 2013
Posts
4
Rep Power
0

#### Rotation angle using 3 points

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).

Code:
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.

[IMG][/IMG]

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.

Thanks!!!
2. This post shows a misunderstanding of the original problem, and was furthermore misunderstood by the thread originator. Disregarding most of this post is a good idea[/edit]

I see two problems. First, atan2 arguments are in order y, x not atan2(x,y). Secondly, where you have
AB=[B[0]-A[0],B[1]-A[1]]
BC=[C[0]-B[0],C[1]-B[1]]
should be
AB=[A[0]-B[0],A[1]-B[1]]
BC=[C[0]-B[0],C[1]-B[1]]
I'll dissuade you from arguing that you've used clever trigonometric identities. When the program doesn't work, simplify. The following analysis shows how I've come to this conclusion.
Code:
2D Points A, B, C look like this when plotted:

B (4,6)
A (-1,4)

C (7,-2)

Vector from B to A named BA would have
length (Ax-Bx) in the \hat{x} direction,
and (Ay-By) in the \hat{y} direction.

The angle is atan2(Ay-By, Ax-Bx).
verify for the sample data.
atan2(-2,-5) is in quadrant 3.  Looks correct.

angle of vector BA: -160   degrees approximately
angle of vector BC:  -70   degrees approximately

You want positive angles.  Add 360 and take modulus 360.
angle BA 200
angle BC 290

The subtended angle by ABC is
(290-200) degrees = 90 degrees
Since it's less than 180 my answer
is 90 degrees counter-clockwise.

If instead we wanted CBA we'd find angle of -90,
adding 360 then mod 360 we get 270
which is greater than 180 so I'd call the angle
from BC to BA as 90 degrees clockwise.

Otherwise,
Code:
for indexA in range(len(coords)-3):
A=a=coords[indexA]
B=coords[indexA+1]
C=coords[indexA+2]
will fix your IndexError. And with a bit of cleanup and cutting before the index error your output looks like
Code:
\$ python p.py
A                       B                       C                       orientation     angle   bearing AB      bearing BC
[283907, 971700]        [284185, 971634]        [284287, 971507]        right             37.87  103.35         141.23
[284185, 971634]        [284287, 971507]        [284275, 971608]        left            -148.00  141.23         353.22
[284287, 971507]        [284275, 971608]        [283919, 971761]        left             -59.96  353.22         293.26
[284275, 971608]        [283919, 971761]        [284311, 971648]        right            172.82  293.26         106.08
[283919, 971761]        [284311, 971648]        [284277, 971637]        right            145.99  106.08         252.07
[284311, 971648]        [284277, 971637]        [284280, 971651]        right            120.02  252.07          12.09
[284277, 971637]        [284280, 971651]        [284174, 971649]        left            -103.17   12.09         268.92
[284280, 971651]        [284174, 971649]        [283909, 971701]        right             12.18  268.91         281.1
[284174, 971649]        [283909, 971701]        [283941, 971700]        right            170.68  281.1           91.79
[283909, 971701]        [283941, 971700]        [284294, 971518]        right             25.48   91.782        117.27
[283941, 971700]        [284294, 971518]        [284288, 971517]        right            143.26  117.27         260.54
[284294, 971518]        [284288, 971517]        [284315, 971539]        right            150.28  260.54          50.83
[284288, 971517]        [284315, 971539]        [284250, 971505]        left            -168.43   50.826        242.39
Which rows did you identify as incorrect?
What should they have read?
Last edited by b49P23TIvg; March 16th, 2013 at 11:53 AM.
3. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Mar 2013
Posts
4
Rep Power
0

#### thanks, but...

Hi! Thanks for your answer.
I will not argue I'm using "clever trigonometric identities", but on the contrary, that they are oretty dumb. If my code is complex it's because I don't remember enough basic trig to make it simpler ;P
I inverted the x-y in the arctan formula because I was calculating the bearing and the angle in that case is the one formed with the Y axis, rather than the X. But then again, I'm not a 100% sure of how arctan works.
Anyways, we both obtained exactly the same results Too bad I didn't realize the table where I showed the difference between my numbers and the ones obtained with GME didn't upload..
Here it is again. The first part is my result, and at the end you can see the ones obtained with that program.

Code:
	A			B			C		orientation	angle		bearing AB	bearing BC	GME angle	GME bearing AB	GME bearing 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

It might be the case that this other program is wrong; if so, can you think of a way of testing what is the real number?

Thanks!!!

Noelia
4. The image you've now inserted is beautiful.

Summary of this result.
• On Cartesian grid I compute same bearings as did noebyus. I used executable Iverson notation, work shown below.
• I computed the discrepancy between bearing and GME bearing, rounded before subtraction to nearest degree.
• I found a strong negative correlation between the discrepancy and the first value of the second pair. Near the end of this post I'll show the meaning. I named the discrepancy ERROR .

I'll describe the problem as I now understand it:
"How do I follow a path?"
Answer:
Code:
orient yourself in start direction,
travel length of first vector.
while unfinished:
turn so many degrees to the {right|left},
travel length of next vector
On a flat Cartesian grid I compute the same bearings as you. How could the GME results be correct? Frankly they aren't making sense to me even on the surface of a sphere way off the equator. If I knew they were spherical coordinates and you explained them I'd try harder. If you read through to the length/bearing table result below you'll see that I checked vector length as a possible contribution to the discrepancy. I couldn't find a correlation.
Code:
   DATA
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

[DIFF =: differences DATA
278  _66
102 _127
_12  101
_356  153
392 _113
_34  _11
3   14
_106   _2
_265   52
32   _1
353 _182
_6   _1
27   22
_65  _34

[COMPLEX =: j./"1 DIFF
278j_66 102j_127 _12j101 _356j153 392j_113 _34j_11 3j14 _106j_2 _265j52 32j_1 353j_182 _6j_1 27j22 _65j_34

[RADIANS =: 12 o. COMPLEX
_0.233095 _0.894138 1.68905 2.73568 _0.280657 _2.82869 1.3597 _3.12273 2.94783 _0.0312398 _0.476034 _2.97644 0.683709 _2.65965

[DEGREES =: degrees RADIANS
_13.3553 _51.2303 96.7757 156.743 _16.0804 _162.072 77.9052 _178.919 168.898 _1.78991 _27.2748 _170.538 39.1737 _152.387

[BEARING =: 90-DEGREES
103.355 141.23 _6.77566 _66.7432 106.08 252.072 12.0948 268.919 _78.8981 91.7899 117.275 260.538 50.8263 242.387

(;:'LENGTH BEARING GME_BEARING'),:,.&.><"1 (round (10 o. COMPLEX) ,: 360 | 360 + BEARING),GME_BEARING
┌──────┬───────┬───────────┐
│LENGTH│BEARING│GME_BEARING│
├──────┼───────┼───────────┤
│286   │103    │103        │
│163   │141    │141        │
│102   │353    │353        │
│387   │293    │284        │
│408   │106    │106        │
│ 36   │252    │252        │
│ 14   │ 12    │ 12        │
│106   │269    │269        │
│270   │281    │138        │
│ 32   │ 92    │  9        │
│397   │117    │ 85        │
│  6   │261    │261        │
│ 35   │ 51    │ 51        │
│ 73   │242    │263        │
└──────┴───────┴───────────┘

differences
2&(-~/\)
round
<.@:(1r2&+)
degrees
57.2957795130823229&*
Correlations:
Code:
   ERROR=:BEARING-GME_BEARING

load'stats'

BEARING corr BEARING  NB. corr gives correlation coefficient
1

LENGTH corr ERROR   NB. vector length doesn't correlate with discrepancy
0.176669

BEARING corr ERROR  NB. bearing doesn't correlate with discrepancy
0.0653035

(}.{.|:DATA) corr ERROR   NB. Possible strong (negative) correlation, remember I'm using rounded values
_0.721059

NB. The first two columns have the strong negative correlation.  BEARING for reference
(}.{.|:DATA) ,.ERROR ,. BEARING
284185   0 103
284287   0 141
284275   0 353
283919   9 293
284311   0 106
284277   0 252
284280   0  12
284174   0 269
283909 143 281
283941  83  92
284294  32 117
284288   0 261
284315   0  51
284250 _21 242
Last edited by b49P23TIvg; March 16th, 2013 at 11:54 AM.
5. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Mar 2013
Posts
4
Rep Power
0
Wow! Thanks for all your work!!!
I'm working on UTM coordinate system, which is not spherical, so that shouldn't be the problem.

It's good to know that the error seems to rely on something related to the first value of the second pair. I'll give another look at the input to see if there's any mistake there.

Noe
6. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Mar 2013
Posts
4
Rep Power
0
Sorry for the delay in making an update on the issue; I had some chaotic weeks and only now can sit to write something down.

So, the issue ended up being my input. I had inadvertently worked with two different files, one of which was in chronological order and one that wasn't. That's why some numbers coincided and some didn't.

In conclusion, both my code, b49P23TIvg's code and GME's code works.

I guess I needed a second pair of eyes to make me realize the code itself wasn't the problem.