Bearing between two points

Question

Ive been using the geopy package , which does a great job, however some of the results i get are inconsistant or come with a relatively large displacement, i suspect that the problem resides with my bearing calculation:

def gb(x,y,center_x,center_y):
dx=x-center_x
dy=y-center_y
if ((dy>=0)and((dx>0)or(dx<0))):
    return math.degrees(math.atan2(dy,dx))
elif (dy<=0)and((dx>0)or (dx<0)):
    return (math.degrees(math.atan2(dy,dx))+360)
else:
    return (math.degrees(math.atan2(dy,dx))+360)%360

I need to calculate the bearing, s.t. center_x and center_y are the pivot. afterwards i use geopy to reverse engineer the gps coordinate:

latlon = VincentyDistance(miles=dist).destination(Point(lat1, lon1), bearing)

Can anyone point me to what might i be doing wrong?

Fiscal

If your problem has been solved by an answer, you should accept the answer by clicking on the big "tick" next to the answer. Otherwise, give feedback so that we can try to help you further. – Germanous

Answer 1

Can anyone point me to what might i be doing wrong?

1. Not showing an example of your "inconsistant or come with a relatively large displacement" results nor your expected results; consequently answerers must rely on guesswork.

2. Not saying what units your input (x, y, etc) is measured in, and how you get the dist used in the destination calculation. I'm assuming (in calculating bearing2 below) that positive x is easting in miles and positive y is northing in miles. It would help greatly if you were to edit your question to fix (1) and (2).

3. A coding style that's not very conducive to making folk want to read it ... have a look through this.

4. In school trigonometry, angles are measured anticlockwise from the X axis (East). In navigation, bearings are measured clockwise from the Y axis (North). See code below. For an example of bearing in use, follow this link, scroll down to the "Destination point given distance and bearing from start point" section, notice that the example is talking about bearings of about 96 or 97 degrees, then click on "view map" and you will notice that the heading is slightly south of east (east being 90 degrees).

Code:

from math import degrees, atan2
    def gb(x, y, center_x, center_y):
        angle = degrees(atan2(y - center_y, x - center_x))
        bearing1 = (angle + 360) % 360
        bearing2 = (90 - angle) % 360
        print "gb: x=%2d y=%2d angle=%6.1f bearing1=%5.1f bearing2=%5.1f" % (x, y, angle, bearing1, bearing2)

    for pt in ((0, 1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1),(-1, 0),(-1,1)):
        gb(pt[0], pt[1], 0, 0)

Output:

gb: x= 0 y= 1 angle=  90.0 bearing1= 90.0 bearing2=  0.0
gb: x= 1 y= 1 angle=  45.0 bearing1= 45.0 bearing2= 45.0
gb: x= 1 y= 0 angle=   0.0 bearing1=  0.0 bearing2= 90.0
gb: x= 1 y=-1 angle= -45.0 bearing1=315.0 bearing2=135.0
gb: x= 0 y=-1 angle= -90.0 bearing1=270.0 bearing2=180.0
gb: x=-1 y=-1 angle=-135.0 bearing1=225.0 bearing2=225.0
gb: x=-1 y= 0 angle= 180.0 bearing1=180.0 bearing2=270.0
gb: x=-1 y= 1 angle= 135.0 bearing1=135.0 bearing2=315.0

Germanous

Hello John, youre right about my coding style, thanks for bringing it to my attention.. as for the code i believe that the implementation youve provided works better, my code for getting a bearing was simply taken form wikipedia, im gonna run my script again and test the results,thanks! – Fiscal

@242Eld: Please supply the URL for the wikipedia article that you mentioned. Also you haven't addressed points (1) and (2) of my answer. – Germanous

the wikipedia article is in hebrew. since im doing a reversed geocoding, an example could not be described with code.. however the problem domain is that after parsing an html page from bing maps, i need to reverse the geocoding process s.t. i could retrieved a desired polygon. the result is good but not satisfying, since there is a displacement of 200-1500 meters from the polygon's original position. im using x y coordinates in pixels, whom i convert to miles using a given scale. the reversed geocoding is made using the geopy package's VincentyDistance – Fiscal

This does not answer the original question but john-machin asked for the wiki site woth the math, this site mentions the same math as 242Eld uses. if anyone should be looking for this. – Patois

Answer 2

I'm not quite sure what you're trying to do in your code,but I do see some oddities that may need to be cleaned up.

1. You test for dy<=0 after you test for dy>=0 in the conditional before. What should your code do if dy==0 and dx==0.
2. Your test ((dy>=0)and((dx>0)or(dx<0))) is equivalent to (dy>=0 and dx!=0), is this what you intended?
3. You are basically doing the same thing in all of your conditionals. Couldn't return math.degrees(math.atan2(dy,dx))+360)%360 work in every scenario? In that case, you wouldn't need to use your if statements anyway.

Goodale

+1 especially for 3, but it may not answer the question. Do you know if e.g. math.atan2 might give relatively inaccurate results near dx==0 (does it include some logic to avoid division-by-nearly-zero issues)? If an if statement is needed, it may be to reflect around a 45 degree line, to get better numeric accuracy - but I would have thought the library would handle that anyway. – Rana

@Steve314, Good point. I looked up the Python2.7 documents and it doesn't gave any warnings concerning an input of 0. Furthermore, I tried both math.atan2(0,30) and math.atan2(30,0) and neither returned an error, which I would expect a built-in library to do if you gave it invalid inputs. – Goodale