4

I would like to find the angle of two lines at their intersection point. In my opinion it is the azimuth I'm looking for, right?

So I'm using Python with shapely for finding the intersection and QGIS for visualisation, but I'm not able to calculate the angle with Python.... So I decided to look how QGIS is doing this

  double QgsPoint::azimuth( const QgsPoint& other ) const
   {
     double dx = other.x() - m_x;
     double dy = other.y() - m_y;
     return ( atan2( dx, dy ) * 180.0 / M_PI );
   }

QGis

But i'm not able to reproduce this result

>>> x2=47.052246
>>> x1=47.052257
>>> y2=8.296885
>>> y1=8.296906
>>> xAns=x1-x2
>>> yAns=y1-y2
>>> np.arctan2(xAns,yAns)*180/np.pi
27.645975364825276

With QGis i'm able to get the angle graphically from Point2 over the intersection too Point1 enter image description here

Furthermore I tried

α = tan-1(cos(lat1) (long2 − long1) / (lat2 − lat1)) 

from ->charlespetzold.com/etc/AvenuesOfManhattan/

without success....

#

Ok after the help of @vagvaf and @vince i wrote the following (not final)

from shapely.geometry import Point
import numpy as np

P1 = Point(47.052257,8.296906)
P2 = Point(47.052246,8.296885)
interP = Point(47.052150,8.296925)

dx = P1.x-interP.x
dy = P1.y-interP.y

dx2 = interP.x-P2.x
dy2 = interP.y-P2.y

azimuth1 = np.arctan2(dx,dy)*180/np.pi
azimuth2 = np.arctan2(dx2,dy2)*180/np.pi

print 'dx ' + str(dx) + '\n' + 'dy ' + str(dy) + '\n' + 'azimuth P1->interP ' + str(azimuth1)
azimuth11 = 360-azimuth1    # dx>0 & dy<0
print '360 - [azimuth1 P1->interP] ' + str(azimuth11)
print 'dx2 ' + str(dx2) + '\n' + 'dy2 ' + str(dy2) + '\n' + 'azimuth interP->P2 ' + str(azimuth2)
azimuth22 = 180-azimuth2    # dx<0 & dy>0
print '180 - [azimuth2 interP->P2] ' + str(azimuth22)
print 'azimuth11 - azimuth22 ' + str(azimuth11-azimuth22)

that results in:

dx 0.000107
dy -1.9e-05
azimuth P1->interP 100.069062699
360 - [azimuth1 P1->interP] 259.930937301
dx2 -9.59999999992e-05
dy2 4.00000000003e-05
azimuth interP->P2 -67.3801350517
180 - [azimuth2 interP->P2] 247.380135052
azimuth11 - azimuth22 = 12.5508022494

that is more or less equal to the qgis result

  • 1
    You can't apply Cartesian trigonometry to angular coordinate values. You must either project into an azmuthal projection or use algorithms developed to solve the "direct problem" of geodesy to determine the bearing. – Vince Jul 5 '16 at 11:02
4

Your script is too complex, simple use the modulo function

def azimuth(point1, point2):
    '''azimuth between 2 shapely points (interval 0 - 360)'''
    angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
    return np.degrees(angle)if angle>0 else np.degrees(angle) + 360

azimuth(interP,P2)
112.61986494834154
azimuth(P2,interP)
292.61986494834156
azimuth(P1,P2)
207.64597536482526
azimuth(P2,P1)
27.645975364825276
azimuth(P1,P2)
207.64597536482526
azimuth(P2,P1)
27.645975364825276
# result
azimuth(P2,interP) - azimuth(P1,interP)
12.550802249443507
# or
azimuth(interP,P2) - azimuth(interP,P1)
12.550802249443507

But you use here the Euclidean distance and if your points are geodetic (angles, longitude, latitude) this distance is meaningless (see inaccurate distance measurements in Python).

With the pygc module (Vincenty's formulae)

print P1.distance(P2) # shapely distance in "angular unit"
2.37065391828006e-05

from pygc import * # distance in meters
print great_distance(start_latitude=P1.y, start_longitude=P1.x,    end_latitude=P2.y,end_longitude=P2.x)
{'reverse_azimuth': 27.553194547658386, 'distance': array(2.619663014755096), 'azimuth': 207.55319613451832}
P2az = great_distance(start_latitude=P2.y, start_longitude=P2.x,    end_latitude=interP.y,end_longitude=interP.x)['reverse_azimuth']
P1az = great_distance(start_latitude=P1.y, start_longitude=P1.x,    end_latitude=interP.y,end_longitude=interP.x)['reverse_azimuth']
print P2az -P1az
12.592170938638404

And as Charles Petzold concludes (How Far from True North are the Avenues of Manhattan? ), the angular difference is minimal

1

In order to calculate the angle P1-InterPoint-P2:

Calculate the difference of the azimuths InterPoint-P1 and Interpoint-P2.

Mind you that you have to investigate the azimuth according to the quadrant that it lies, so:

if dx>0 & dy>0 : final azimith=azimuth
if dx<0 & dy>0 : final azimuth=180-azimuth
if dx<0 & dy<0: final azimuth=180+azimuth
if dx>0 & dy<0 final azimuth=360-azimuth
  • Thank you for your answer vagvaf. When i calculate the diff of the azimuth i get the exterior angle, am i right? so i subtracted the result from 180 and get the angle for the interior (the same result gives me qgis).is this correct or more by accident? – Hendrik Jul 5 '16 at 10:48
  • a simple subtraction would give you the result. Take into consideration that, when calculating azimuths you have to investigate according to the quadrant that your azimuth lies. – vagvaf Jul 5 '16 at 11:06
  • It's not stated, but the coordinate values are geodetic (either that, or submillimeter resolution). This answer is more like a comment, in that it doesn't explain how to do what it suggests. – Vince Jul 5 '16 at 11:13
  • @vagvaf the second and third if statement are equal. if dx<0 & dy>0 : final azimuth=180-azimuth, if dx>0 & dy<0 : final azimuth=180+azimuth – Hendrik Jul 5 '16 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.