--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/geometry.py Wed Jul 29 23:45:53 2020 +0300
@@ -0,0 +1,124 @@
+from math import sqrt, hypot, radians, sin, cos, atan2
+
+class Location:
+ def __init__(self, latitude, longitude):
+ self.latitude, self.longitude = latitude, longitude
+ def distance(self, other):
+ # https://stackoverflow.com/a/365853
+ a = sin(radians(self.latitude - other.latitude) / 2) ** 2
+ a += sin(radians(self.longitude - other.longitude) / 2) ** 2 * cos(radians(self.latitude)) * cos(radians(other.latitude))
+ return 6371 * 2 * atan2(sqrt(a), sqrt(1 - a))
+ def __repr__(self):
+ return '%s(%r, %r)' % (type(self).__name__, self.latitude, self.longitude)
+ def __str__(self):
+ return '%.5f, %.5f' % (self.latitude, self.longitude)
+ @property
+ def x(self):
+ return self.latitude
+ @property
+ def y(self):
+ return self.longitude
+ @property
+ def link_to_map(self):
+ return 'http://www.openstreetmap.org/#map=19/%f/%f' % (self.latitude, self.longitude)
+
+class Ring:
+ def __init__(self, container):
+ self.container = container
+ def __getitem__(self, i):
+ while i < 0:
+ i += len(self.container)
+ while i >= len(self.container):
+ i -= len(self.container)
+ return self.container[i]
+ def __iter__(self):
+ return iter(self.container)
+ def __len__(self):
+ return len(self.container)
+
+class Polygon:
+ def __init__(self, *points):
+ self.points = points
+ def __repr__(self):
+ return '%s(%s)' % (type(self).__name__, ', '.join(map(repr, self.points)))
+ def area(self):
+ ring = Ring(self.points)
+ return sum(
+ ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y
+ for i in range(len(ring))
+ ) / 2
+ def circumference(self):
+ ring = Ring(self.points)
+ return sum(
+ sqrt((ring[i + 1].x - ring[i].x)**2 + (ring[i + 1].y - ring[i].y)**2)
+ for i in range(len(ring))
+ )
+ def centroid(self):
+ ring = Ring(self.points)
+ x = sum(
+ (ring[i].x + ring[i + 1].x) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y)
+ for i in range(len(ring))
+ ) / 6 / self.area()
+ y = sum(
+ (ring[i].y + ring[i + 1].y) * (ring[i].x * ring[i + 1].y - ring[i + 1].x * ring[i].y)
+ for i in range(len(ring))
+ ) / 6 / self.area()
+ return self.point_type()(x, y)
+ def point_type(self):
+ if len(self.points):
+ return type(self.points[0])
+ else:
+ return Point
+ def segments(self):
+ ring = Ring(self.points)
+ for i in range(len(ring)):
+ yield LineSegment(ring[i], ring[i + 1])
+ def contains_point(self, point):
+ outer_point = self.point_type()(
+ min(point.x for point in self.points) - 1,
+ min(point.y for point in self.points) - 1
+ )
+ outer_segment = LineSegment(point, outer_point)
+ intersections = 0
+ for segment in self.segments():
+ if segment.intersection(outer_segment) is not None:
+ intersections += 1
+ return bool(intersections & 1)
+
+class LineSegment:
+ def __init__(self, p1, p2):
+ self.p1, self.p2 = p1, p2
+ def __repr__(self):
+ return 'LineSegment(%r, %r)' % (self.p1, self.p2)
+ def length(self):
+ return hypot(self.p1.x - self.p2.x, self.p1.y - self.p2.y)
+ def intersection(self, other):
+ point_type = type(self.p1)
+ x = (self.p1.x, self.p2.x, other.p1.x, other.p2.x)
+ y = (self.p1.y, self.p2.y, other.p1.y, other.p2.y)
+ try:
+ denominator = (x[0] - x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] - x[3])
+ Px = ((x[0] * y[1] - y[0] * x[1]) * (x[2] - x[3]) - (x[0] - x[1]) * (x[2] * y[3] - y[2] * x[3])) / denominator
+ Py = ((x[0] * y[1] - y[0] * x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] * y[3] - y[2] * x[3])) / denominator
+ distance = lambda n: hypot(Px - x[n], Py - y[n])
+ if max(distance(0), distance(1)) <= self.length() and max(distance(2), distance(3)) <= other.length():
+ return point_type(Px, Py)
+ else:
+ return None
+ except ZeroDivisionError:
+ return None
+
+class Point:
+ def __init__(self, x, y):
+ self.x, self.y = x, y
+ def __repr__(self):
+ return 'Point(%r, %r)' % (self.x, self.y)
+
+A = Polygon(
+ Point(2,3),
+ Point(1,1),
+ Point(4,0),
+ Point(6,2),
+ Point(4,4))
+L1 = LineSegment(Point(1, 1), Point(-1, 5))
+L2 = LineSegment(Point(1, 5), Point(5, 1))
