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.