[COINS] angles between segments in the same quadrant of the Cartesian plane

[fnattino]

Describe the problem

Hi there, thanks a lot for this fantastic library!

I am playing with momepy's COINS implementation and I have identified some strange behaviour when creating strokes over edges that connect with some sharp angles. I have narrowed down the issue to what I believe being wrong angles being calculated for segments lying in the same quadrant of the Cartesian plane (see examples below).

If I understand _angle_between_two_lines correctly, one checks whether two segments lie on the same side of the horizontal plane (e.g. here). If this is found to be the case, the result seems to assume that the two segments lie on opposite sides of the vertical plane (angle = 180 - (abs(l1orien) + abs(l2orien))). The case with the segments lying on the same side of the vertical plane (angle = abs(l1orien - l2orien)) seems not to be covered. Or am I missing something here?

Happy to contribute with a PR if you confirm this is unwanted behaviour!

Steps to reproduce

For the following pair of segments, I get these angles:

import momepy.coins

l1 = [[0., 0.], [0.1, 0.2]]
l2 = [[0., 0.], [0.2, 0.1]]
angle_l12 = momepy.coins._angle_between_two_lines(l1, l2)
print(angle_l12)
# 90.0

l3 = [[-0.4, -0.1], [0., 0.]]
l4 = [[-0.6, -0.2], [0., 0.]]
angle_l34 = momepy.coins._angle_between_two_lines(l3, l4)
print(angle_l34)
# 147.529

while I would have expect instead something like:

import numpy as np

def normalize(x):
    return np.array(x) / np.linalg.norm(x)

def angle_between(l1, l2):
    v1 = np.array(l1)
    v1 -= v1[0]
    v2 = np.array(l2)
    v2 -= v2[0]
    return np.degrees(np.arccos(np.dot(
        normalize(v1[1]), 
        normalize(v2[1])
    )))

angle_l12 = angle_between(l1, l2)
print(angle_l12)
# 36.86989764584403

angle_l34 = angle_between(l3, l4)
print(angle_l34)
# 4.398705354995591

Visualization of the segments:

Versions of your packages

momepy=0.6.0

Your operating system

macOS - x86-64

Additional context

No response

[martinfleis]

Hey, we've actually noticed the same last week and @anastassiavybornova might be looking into that as well. If you have a fix, PR would be very welcome!

[anastassiavybornova]

yes, looking into it as we speak, i also suspect there is a bug in the internal _points_set_angle function, happy to look further into what you have to share @fnattino on this, too! (in particular, i think the cases D and E from Table 1 here: http://journals.sagepub.com/doi/10.1177/2399808320967680 are not correctly identified by the mentioned func)

[anastassiavybornova]

...and in fact yes, the case you're describing (both lines in the same quadrant) is not even accounted for as far as i can tell

[fnattino]

Hi @anastassiavybornova @martinfleis, great to hear you also noticed this! This crude fix here seems to solve the problem: https://github.com/pysal/momepy/compare/main...fnattino:momepy:fix-angle-between-lines?expand=1 However, I wonder whether one could do this more cleanly by getting _points_set_angle to work for all cases..

[martinfleis]

Maybe tangential but we also have a function to measure angles here https://github.com/pysal/momepy/blob/main/momepy/functional/_dimension.py

[anastassiavybornova]

Maybe tangential but we also have a function to measure angles here https://github.com/pysal/momepy/blob/main/momepy/functional/_dimension.py

@martinfleis just to be sure, you mean _get_angle_njit() here?

[martinfleis]

Yes. I was on the phone so just wanted to drop that in without having time to look at it. It seems to be the same formula used in existing coins and in #630 as well, so it was a pointless suggestion :).