2

the function below tries to calculate the difference between 2 angles. I use the mathematical reasoning is as follows :

1) Obtain the 2 angles :

angle1 = ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1))
angle2 = ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2))

2) Obtain vectors (of lenght == 1) from 2 angles :

V1 = (cos(angle1), sin(angle1))
V2 = (cos(angle2), sin(angle2))

3) Compute the dot product :

V1 * V2 = V1[x] * V2[x] + V1[y] * V2[y]
=
cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1)

Since the lenght of the vectors are both 1, it follows that :

the angle between the two vectors = ArcCos(dotProduct)/(length(v1)*legnth(v2)

= ArcCos(dotProduct)

The problem in the function below is that the dot product is yielding values greater than 1, which should not be possible if the math reasoning is correct, and that is my question :

is the reasoning correct ?

CREATE OR REPLACE FUNCTION angleDiff(l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
BEGIN           
  SELECT ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1)) into angle1;
  SELECT ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2)) into angle2;  

  RETURN acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1));
END;
$$  LANGUAGE plpgsql
  • There appears to be a typographical error in (3) and in the code: one of the sin(angle1) should be sin(angle2). Is that the source of the problem or does fixing that still produce an error? Incidentally, most people would compute the difference between two angles as angle1-angle2. What is the reason you are using this roundabout calculation? – whuber Nov 26 '13 at 17:48
  • the angles 359 and 0 are only one degree appart, but 359 - 0 = 359. ArcCos(dotProduct) will give 1, which is what I want. – Max L. Nov 26 '13 at 17:57
  • 1
    If you would make your requirements more precise, readers could suggest better solutions. For instance, if you want an unsigned angle between 0 and 180 degrees, this can be accomplished with simple arithmetical operations on the angles themselves without the detour through trig and inverse trig functions. For instance, 1 = abs(mod(359-0,360) - 180). – whuber Nov 26 '13 at 20:08
  • Ok, so I take note of your comment regarding the precision of the requirement. As for your solution, it is much more elegant, if you post is as an answer I will vote for it. – Max L. Nov 26 '13 at 20:48
  • My first thought was to use my answer below. Funny thing is, however, i did try using the mod() function, and the % operator, but either way, i couldn't get passed the "integer vs real" error! Some form of casting was needed, but i ended up with my first idea. – Martin F Nov 28 '13 at 23:23
5

This merges your logic and whuber's logic, except the return is signed:

CREATE OR REPLACE FUNCTION angleDiff (l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
DECLARE diff FLOAT;
BEGIN           
    SELECT ST_Azimuth (ST_StartPoint(l1), ST_EndPoint(l1)) INTO angle1;
    SELECT ST_Azimuth (ST_StartPoint(l2), ST_EndPoint(l2)) INTO angle2;  
    SELECT degrees (angle2 - angle1) INTO diff;
    CASE
    WHEN diff >   180 THEN RETURN diff - 360;
    WHEN diff <= -180 THEN RETURN diff + 360;
    ELSE                   RETURN diff;
    END CASE;
END;
$$  LANGUAGE plpgsql

The return value is clockwise-positive angle [0 to 180] from line1 to line2. If line2 is "clockwise before" line1, the return value is negative (0 to -180). Ignore the sign if you don't care about order.

  • the result of st_azimuth is in radian, so it must be converted to degrees (there is the degrees() function for that), otherwise it looks good. – Max L. Nov 29 '13 at 1:03
  • Thanks -- degrees() has now been added. I had only tested a version of angleDiff() that took in two angles, in degrees, and forgot that st_azimuth returns radians. – Martin F Nov 29 '13 at 1:14
  • ok, so it deserves a vote and and "accepted anwer" ;-) – Max L. Nov 29 '13 at 14:14
0

Ok, there were 2 problems, but the math reasoning is correct :

1) the dot product was

acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1));
                                                   ^^^^^^^

it should be :

-- should be : 
acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle2));
                                                   ^^^^^^^

and also, the dot product can go slightly outside the range of [-1, 1], so it needs to be checked.

The corrected versions :

CREATE OR REPLACE FUNCTION angleDiff(l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
DECLARE dotProduct FLOAT;
BEGIN           
  SELECT ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1)) into angle1;
  SELECT ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2)) into angle2;  

  select (cos(angle1) * cos(angle2) + sin(angle1) * sin(angle2)) into dotProduct;  

  if dotProduct > 1 then return 0;
  elseif dotProduct < -1 then return 180;
  else return degrees(acos(dotProduct));
  end if;

END;
$$  LANGUAGE plpgsql

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.