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 );
}
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
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
Best Answer
Your script is too complex, simply use the modulo function:
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)
And as Charles Petzold concludes (How Far from True North are the Avenues of Manhattan? ), the angular difference is minimal