Challenge: fix angle calculation to identify un-squared buildings

[emi420]

Do you have good math and/or GIS skills?

We need your help :)

Hi everyone! as part of the HOT efforts to improve the quality of the mapping, we're working on tools for validation.

Challenge: un-squared buildings

One of the things that we're doing is to identify buildings that have not angles squared.

Goals

1. Calculate the correct internal angles of a polygon.
2. Use them to detect if a polygon has angles outside a threshold (ex: 85º to 95º)

Problem

The function being used for calculating each internal angles is the following.

But the calculations are not working as expected.

'''
  Calculate the angle between 3 points
'''
def angle(x1,y1,x2,y2,x3,y3):
    ba = [x1 - x2, y1 - y2]
    bc = [x3 - x2, y3 - y2]
    
    dot_p = ba[0] * bc[0] + ba[1] * bc[1]
        
    cosine_angle = dot_p / (
        ((ba[0] * ba[0] + ba[1] * ba[1]) ** 0.5) * 
        ((bc[0] * bc[0] + bc[1] * bc[1]) ** 0.5)
    )

    angle = math.acos(cosine_angle)
    return round((angle * 180 / math.pi), 1)

For example, if we take this polygon and calculate the angle of each of its 10 corners. All of them are ~90 degrees.

1. [0, 1, 2] 90.4
2. [1, 2, 3] 89.6
3. [2, 3, 4] 90.4
4. [3, 4, 5] 89.6
5. [4, 5, 6] 89.6
6. [5, 6, 7] 90.4
7. [6, 7, 8] 90.4
8. [7, 8, 9] 89.6
9. [8, 9, 0] 89.6
10. [9, 0, 1] 89.6

Is squared? True

But if we take the same polygon and rotated it a little bit, the results change and none of the calculated angles are ~90 degrees.

1. [0, 1, 2] 101.9
2. [1, 2, 3] 78.1
3. [2, 3, 4] 102.0
4. [3, 4, 5] 78.2
5. [4, 5, 6] 78.2
6. [5, 6, 7] 102.0
7. [6, 7, 8] 101.9
8. [7, 8, 9] 78.1
9. [8, 9, 0] 78.1
10. [9, 0, 1] 78.1

Is squared? False

Solution

We invite contributors to find a solution to this problem.

To help testing and showing your solution, there's a Jupyter Notebook available.

You can check the Gist on my account, that you can also download and use on your local Jupyter server:

https://gist.github.com/emi420/32f33dc16e6e1789cd8077306fac53a5

And there's also a copy of it on Google Colab:

https://colab.research.google.com/drive/1ogPJOIiqhE3yi9mexKaCVCR89DQqK6YA#scrollTo=YWzOVVyb4Aaa

Can you help? please post your solution here. It can be a link to your own code, a copy of the notebook or just a comment explaining what's wrong and how can be fixed.

Thanks!

[emi420]

It looks like the problem was with... map projection :)

Adding a function to convert each coordinate from Epsg4326 to Epsg3857 would solve the issue.

This is the C++ version of that function:

void Geo::epsg4326toEpsg3857(double& x, double& y) {
  x = (x * 20037508.34) / 180;
  y = (std::log(std::tan(((90 + y) * M_PI) / 360))) /
      (M_PI / 180);
  y = (y * 20037508.34) / 180;
}

Thanks @caprieldeluca for your help!

[emi420]

Issue was solved, thanks for your help!