#1
  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. #2
  3. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,841
    Rep Power
    480
    [edit]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.
    [code]Code tags[/code] are essential for python code and Makefiles!
  4. #3
  5. 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
  6. #4
  7. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,841
    Rep Power
    480
    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.
    [code]Code tags[/code] are essential for python code and Makefiles!
  8. #5
  9. 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
  10. #6
  11. 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.

IMN logo majestic logo threadwatch logo seochat tools logo