Search code examples
pythoncoordinatesangleqgis

How come my full circle isn't 360 degrees?


I am stuck with a confusing problem. Here's a little background:

I'm working on qgis/python with coordinate points in Lambert93: one central point (my dict key) and several other points gravitating around it. To simplify, the code I've put down just one example:

import numpy as np
import math

dict = {(355385,6.68906e+06): [(355277,6.68901e+06), (355501,6.68912e+06), (355364,6.6891e+06), (355277,6.68901e+06)]}

for key, values in dict.iteritems():
    anglist =[]
    print key
    i=0
    j=1
    for sides in values[:-1]:

        A = np.array(dict[key][i])
        B = np.array(key)
        C = np.array(dict[key][j])

        BA = A - B
        BC = C - B

        cosine_angle = np.vdot(BA, BC) / (np.linalg.norm(BA) * np.linalg.norm(BC))
        angle = (np.degrees(np.arccos(cosine_angle)))

        i+=1
        j+=1 

        anglist.append(angle)
        s = sum(anglist)
    dict[key]= [values, anglist, s] 
print dict

results are :

{(355385, 6689060.0): [[(355277, 6689010.0), (355501, 6689120.0), (355364, 6689100.0), (355277, 6689010.0)], [177.4925133253854, 90.349597027985112, 87.142916297400205], 354.98502665077069]}

As you can see, sum = 354. I have a large set of data and sometimes I get the correct 360, but for the most part I don't. Yet in all logic, by turning around a single point and ending the calculation where it started, the only result i should get is 360.

I have tried a second way just to see if the cosine-angle and angle weren't the problem :

from math import sqrt
from math import acos
import numpy
def angle(a, b, c):

    # Create vectors from points
    ba = [ aa-bb for aa,bb in zip(a,b) ]
    bc = [ cc-bb for cc,bb in zip(c,b) ]

# Normalize vector
    nba = sqrt ( sum ( (x**2.0 for x in ba) ) )
    ba = [ x/nba for x in ba ]

    nbc = sqrt ( sum ( (x**2.0 for x in bc) ) )
    bc = [ x/nbc for x in bc ]

# Calculate scalar from normalized vectors
    scale = sum ( (aa*bb for aa,bb in zip(ba,bc)) )

# calculate the angle in radian
    angle = numpy.degrees(acos(scale))
    return angle

print angle((355277,6.68901e+06),(355385,6.68906e+06), (355501,6.68912e+06))
print angle((355501,6.68912e+06),(355385,6.68906e+06), (355364,6.6891e+06))
print angle((355364,6.6891e+06),(355385,6.68906e+06), (355277,6.68901e+06))

But the results are still:

177.492513325
90.349597028
87.1429162974

So I think we can cross the math out of the problem... So one possibility is a problem with how qgis (or python?) manages the coordinates. How can I go around this?

I should say, the codes are largely the same as here, here and here


Solution

  • Thanks to Hellmar's indication, here's the working code. I did two things : gave all the central points (0,0) in my plane, substracted the distance between the central point and (0,0) from all the other related points. I was then able to use atan2 because all the angles I needed to calculate were at (0,0). Atan2 is very practical in this context because it only needs two points to calculate the angle : the angle measured is always the one at 0,0. This probably means I don't have to set my central point as 0,0 as it just get kicked out of the equation. Here's my code. Any further advice is much welcome.

    import numpy as np
    dictionary = {(355385,6.68906e+06): [(355277,6.68901e+06), (355501,6.68912e+06), (355364,6.6891e+06), (355277,6.68901e+06)], (355364,6.6891e+06): [(355261,6.68905e+06), (355385,6.68906e+06), (355481,6.68916e+06), (355340,6.68915e+06), (355261,6.68905e+06)], (355340,6.68915e+06): [(355238,6.68909e+06), (355364,6.6891e+06), (355452,6.68921e+06), (355238,6.68909e+06)]}
    def angle_between(p1, p2):
        ang1 = np.arctan2(*p1[::-1])
        ang2 = np.arctan2(*p2[::-1])
        return np.rad2deg((ang1 - ang2) % (2 * np.pi))
    
    zlist=[]
    newdict={}
    for key, values in dictionary.iteritems():
        xlist =[]
        ylist = []
        i=0
    #print key
    #print key[0]
        for sides in values:
    
            A = dictionary[key][i][0]
            B = key[0]
            C = dictionary[key][i][1]
            D= key[1]
            E = (0.0, 0.0)
    
            o1 = A-B
            o2  = C-D
    
            xlist.append(o1)
            ylist.append(o2)
            ziplist = zip(xlist, ylist)
            i+=1
        zlist.append(ziplist)  
    
    #print dict[key][i][0]
    
    newdict=zlist
    print newdict
    angledict = []
    
    for p in newdict:
        i=0
        j=1
        print p
        for q in p[:-1]:
            A=p[i]
            B=p[j]
            print "A=",A
            print "B=", B
    
            angledict.append(angle_between(A,B))
            i+=1
            j+=1
    
    print angledict
    

    I took the angle_between function from here