0

I would like to check if two polygons built from coordinates touch or overlap. The following example should touch (they came from a simulation), however, the code doesn't detect they touch, hence, crash. Any clue?

import geopy
import geopy.distance
import math
from shapely import geometry
def bearing_calc(a_lat, a_lon, b_lat, b_lon):  # a previous position b current position
    # print(a)
    # a.type
    rlat1 = math.radians(a_lat)
    # print(a_lat)
    rlon1 = math.radians(a_lon)
    rlat2 = math.radians(b_lat)
    rlon2 = math.radians(b_lon)
    dlon = math.radians(b_lon - a_lon)
    b = math.atan2(math.sin(dlon) * math.cos(rlat2),
                   math.cos(rlat1) * math.sin(rlat2) - math.sin(rlat1) * math.cos(rlat2) * math.cos(
                       dlon))  # bearing calc
    bd = math.degrees(b)
    br, bn = divmod(bd + 360, 360)  # the bearing remainder and final bearing
    return bn

# To get a rotated rectangle at a bearing, you need to get the points of the the recatangle at that bearing
def get_rotated_points(coordinates, bearing, width, length):
    start = geopy.Point(coordinates)
    width = width / 1000
    length = length / 1000
    rectlength = geopy.distance.distance(kilometers=length)
    rectwidth = geopy.distance.distance(kilometers=width)
    halfwidth = geopy.distance.distance(kilometers=width / 2)
    halflength = geopy.distance.distance(kilometers=length / 2)

    pointAB = halflength.destination(point=start, bearing=bearing)
    pointA = halfwidth.destination(point=pointAB, bearing=0 - bearing)
    pointB = rectwidth.destination(point=pointA, bearing=180 - bearing)
    pointC = rectlength.destination(point=pointB, bearing=bearing - 180)
    pointD = rectwidth.destination(point=pointC, bearing=0 - bearing)

    points = []
    for point in [pointA, pointB, pointC, pointD]:
        coords = (point.latitude, point.longitude)
        points.append(coords)

    return points

v1_id = 176
v1_coord_5 = [41.39129840757358, 2.1658667401858738]

v1_coord_4 = [41.391335226146495, 2.1658363842310404]
v1_length = 5 #meters
v1_width = 1.8 #meters
v1_bearing = bearing_calc(v1_coord_4[0], v1_coord_4[1],
                          v1_coord_5[0], v1_coord_5[1])
v1_points = get_rotated_points(tuple(v1_coord_5), v1_bearing,
                               v1_width, v1_length)
polygon1 = geometry.Polygon(v1_points)

v2_id = 210

v2_coord_5 = [41.39134056977012, 2.1658847515529804]

v2_coord_4 = [41.39127485461452, 2.165851515431892]

v2_length = 5 #meters
v2_width = 1.8 #meters
v2_bearing = bearing_calc(v2_coord_4[0], v2_coord_4[1],
                          v2_coord_5[0], v2_coord_5[1])
v2_points = get_rotated_points(tuple(v2_coord_5), v2_bearing,
                               v2_width, v2_length)
polygon2 = geometry.Polygon(v2_points)

if polygon1.intersection(polygon2).area > 0.0:
  print("COLISION")
else:
  print("NO COLISION")

I've also tried the function touches

polygon1.touches(polygon2)

Unfortunately, also, false. I don't know if there is a function to check if the two borders touch even in one points, minimally.

If I plot the two polygons, they do not collide, as shown above, however, they should. Is it because the code is not considering the earth surface curvature? How can I adapt the code to consider that factor.

enter image description here

GoT GOt
  • 49
  • 1
  • 9

0 Answers0