#!/usr/bin/env python # -*- coding: utf-8 -*- # Copyright (c) 2010. # SMHI, # Folkborgsvägen 1, # Norrköping, # Sweden # Author(s): # Martin Raspaud # This file is part of mpop. # mpop is free software: you can redistribute it and/or modify it under the # terms of the GNU General Public License as published by the Free Software # Foundation, either version 3 of the License, or (at your option) any later # version. # mpop is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with # mpop. If not, see . """Spherical geometry module. """ import math import numpy as np EPSILON = 0.0000001 class Coordinate(object): """Point on earth in terms of lat and lon. """ lat = None lon = None x__ = None y__ = None z__ = None def __init__(self, lat=None, lon=None, x__=None, y__=None, z__=None, R=1): self.R = R if lat is not None and lon is not None: self.lat = math.radians(lat) self.lon = math.radians(lon) self._update_cart() else: self.x__ = x__ self.y__ = y__ self.z__ = z__ self._update_lonlat() def _update_cart(self): """Convert lon/lat to cartesian coordinates. """ self.x__ = math.cos(self.lat) * math.cos(self.lon) self.y__ = math.cos(self.lat) * math.sin(self.lon) self.z__ = math.sin(self.lat) def _update_lonlat(self): """Convert cartesian to lon/lat. """ self.lat = math.degrees(math.asin(self.z__ / self.R)) self.lon = math.degrees(math.atan2(self.y__, self.x__)) def __ne__(self, other): if(abs(self.lat - other.lat) < EPSILON and abs(self.lon - other.lon) < EPSILON): return 0 else: return 1 def __eq__(self, other): return not self.__ne__(other) def __str__(self): return str((math.degrees(self.lat), math.degrees(self.lon))) def __repr__(self): return str((math.degrees(self.lat), math.degrees(self.lon))) def cross2cart(self, point): """Compute the cross product, and convert to cartesian coordinates (assuming radius 1). """ lat1 = self.lat lon1 = self.lon lat2 = point.lat lon2 = point.lon res = Coordinate( x__=(math.sin(lat1 - lat2) * math.sin((lon1 + lon2) / 2) * math.cos((lon1 - lon2) / 2) - math.sin(lat1 + lat2) * math.cos((lon1 + lon2) / 2) * math.sin((lon1 - lon2) / 2)), y__=(math.sin(lat1 - lat2) * math.cos((lon1 + lon2) / 2) * math.cos((lon1 - lon2) / 2) + math.sin(lat1 + lat2) * math.sin((lon1 + lon2) / 2) * math.sin((lon1 - lon2) / 2)), z__=(math.cos(lat1) * math.cos(lat2) * math.sin(lon1 - lon2))) return res def distance(self, point): """Vincenty formula. """ dlambda = self.lon - point.lon num = ((math.cos(point.lat) * math.sin(dlambda)) ** 2 + (math.cos(self.lat) * math.sin(point.lat) - math.sin(self.lat) * math.cos(point.lat) * math.cos(dlambda)) ** 2) den = (math.sin(self.lat) * math.sin(point.lat) + math.cos(self.lat) * math.cos(point.lat) * math.cos(dlambda)) return math.atan2(math.sqrt(num), den) def norm(self): """Return the norm of the vector. """ return math.sqrt(self.x__ ** 2 + self.y__ ** 2 + self.z__ ** 2) def normalize(self): """normalize the vector. """ norm = self.norm() self.x__ /= norm self.y__ /= norm self.z__ /= norm return self def cross(self, point): """cross product with another vector. """ x__ = self.y__ * point.z__ - self.z__ * point.y__ y__ = self.z__ * point.x__ - self.x__ * point.z__ z__ = self.x__ * point.y__ - self.y__ * point.x__ return Coordinate(x__=x__, y__=y__, z__=z__) def dot(self, point): """dot product with another vector. """ return (self.x__ * point.x__ + self.y__ * point.y__ + self.z__ * point.z__) class Arc(object): """An arc of the great circle between two points. """ start = None end = None def __init__(self, start, end): self.start, self.end = start, end def center_angle(self): """Angle of an arc at the center of the sphere. """ val = (math.cos(self.start.lat - self.end.lat) + math.cos(self.start.lon - self.end.lon) - 1) if val > 1: val = 1 elif val < -1: val = -1 return math.acos(val) def __eq__(self, other): if(self.start == other.start and self.end == other.end): return 1 return 0 def __ne__(self, other): return not self.__eq__(other) def __str__(self): return str((str(self.start), str(self.end))) def angle(self, other_arc): """Oriented angle between two arcs. """ if self.start == other_arc.start: a__ = self.start b__ = self.end c__ = other_arc.end elif self.start == other_arc.end: a__ = self.start b__ = self.end c__ = other_arc.start elif self.end == other_arc.end: a__ = self.end b__ = self.start c__ = other_arc.start elif self.end == other_arc.start: a__ = self.end b__ = self.start c__ = other_arc.end else: raise ValueError("No common point in angle computation.") ua_ = a__.cross(b__) ub_ = a__.cross(c__) val = ua_.dot(ub_) / (ua_.norm() * ub_.norm()) if abs(val - 1) < EPSILON: angle = 0 elif abs(val + 1) < EPSILON: angle = math.pi else: angle = math.acos(val) n__ = ua_.normalize() if n__.dot(c__) > 0: return -angle else: return angle def intersections(self, other_arc): """Gives the two intersections of the greats circles defined by the current arc and *other_arc*. """ if self.end.lon - self.start.lon > math.pi: self.end.lon -= 2 * math.pi if other_arc.end.lon - other_arc.start.lon > math.pi: other_arc.end.lon -= 2 * math.pi if self.end.lon - self.start.lon < -math.pi: self.end.lon += 2 * math.pi if other_arc.end.lon - other_arc.start.lon < -math.pi: other_arc.end.lon += 2 * math.pi ea_ = self.start.cross2cart(self.end).normalize() eb_ = other_arc.start.cross2cart(other_arc.end).normalize() cross = ea_.cross(eb_) lat = math.atan2(cross.z__, math.sqrt(cross.x__ ** 2 + cross.y__ ** 2)) lon = math.atan2(-cross.y__, cross.x__) return (Coordinate(math.degrees(lat), math.degrees(lon)), Coordinate(math.degrees(-lat), math.degrees(modpi(lon + math.pi)))) def intersects(self, other_arc): """Says if two arcs defined by the current arc and the *other_arc* intersect. An arc is defined as the shortest tracks between two points. """ for i in self.intersections(other_arc): a__ = self.start b__ = self.end c__ = other_arc.start d__ = other_arc.end ab_ = a__.distance(b__) cd_ = c__.distance(d__) if(abs(a__.distance(i) + b__.distance(i) - ab_) < EPSILON and abs(c__.distance(i) + d__.distance(i) - cd_) < EPSILON): return True return False def intersection(self, other_arc): """Says where, if two arcs defined by the current arc and the *other_arc* intersect. An arc is defined as the shortest tracks between two points. """ for i in self.intersections(other_arc): a__ = self.start b__ = self.end c__ = other_arc.start d__ = other_arc.end ab_ = a__.distance(b__) cd_ = c__.distance(d__) if(abs(a__.distance(i) + b__.distance(i) - ab_) < EPSILON and abs(c__.distance(i) + d__.distance(i) - cd_) < EPSILON): return i return None def modpi(val): """Puts *val* between -pi and pi. """ return (val + math.pi) % (2 * math.pi) - math.pi def modpi2(val): """Puts *val* between 0 and 2pi. """ return val % (2 * math.pi) def point_inside(point, corners): """Is a point inside the 4 corners ? This uses great circle arcs as area boundaries. """ arc1 = Arc(corners[0], corners[1]) arc2 = Arc(corners[1], corners[2]) arc3 = Arc(corners[2], corners[3]) arc4 = Arc(corners[3], corners[0]) arc5 = Arc(corners[1], point) arc6 = Arc(corners[3], point) angle1 = modpi(arc1.angle(arc2)) angle1bis = modpi(arc1.angle(arc5)) angle2 = modpi(arc3.angle(arc4)) angle2bis = modpi(arc3.angle(arc6)) return (np.sign(angle1) == np.sign(angle1bis) and abs(angle1) > abs(angle1bis) and np.sign(angle2) == np.sign(angle2bis) and abs(angle2) > abs(angle2bis)) def overlaps(area_corners, segment_corners): """Are two areas overlapping ? This uses great circle arcs as area boundaries. """ for i in area_corners: if point_inside(i, segment_corners): return True for i in segment_corners: if point_inside(i, area_corners): return True area_arc1 = Arc(area_corners[0], area_corners[1]) area_arc2 = Arc(area_corners[1], area_corners[2]) area_arc3 = Arc(area_corners[2], area_corners[3]) area_arc4 = Arc(area_corners[3], area_corners[0]) segment_arc1 = Arc(segment_corners[0], segment_corners[1]) segment_arc2 = Arc(segment_corners[1], segment_corners[2]) segment_arc3 = Arc(segment_corners[2], segment_corners[3]) segment_arc4 = Arc(segment_corners[3], segment_corners[0]) for i in (area_arc1, area_arc2, area_arc3, area_arc4): for j in (segment_arc1, segment_arc2, segment_arc3, segment_arc4): if i.intersects(j): return True return False def get_intersections(b__, boundaries): """Get the intersections of *b__* with *boundaries*. Returns both the intersection coordinates and the concerned boundaries. """ intersections = [] bounds = [] for other_b in boundaries: inter = b__.intersection(other_b) if inter is not None: intersections.append(inter) bounds.append(other_b) return intersections, bounds def get_first_intersection(b__, boundaries): """Get the first intersection on *b__* with *boundaries*. """ intersections, bounds = get_intersections(b__, boundaries) del bounds dists = np.array([b__.start.distance(p__) for p__ in intersections]) indices = dists.argsort() if len(intersections) > 0: return intersections[indices[0]] return None def get_next_intersection(p__, b__, boundaries): """Get the next intersection from the intersection of arcs *p__* and *b__* along segment *b__* with *boundaries*. """ new_b = Arc(p__, b__.end) intersections, bounds = get_intersections(new_b, boundaries) dists = np.array([b__.start.distance(p2) for p2 in intersections]) indices = dists.argsort() if len(intersections) > 0 and intersections[indices[0]] != p__: return intersections[indices[0]], bounds[indices[0]] elif len(intersections) > 1: return intersections[indices[1]], bounds[indices[1]] return None, None def polygon(area_corners, segment_corners): """Get the intersection polygon between two areas. """ area_boundaries = [Arc(area_corners[0], area_corners[1]), Arc(area_corners[1], area_corners[2]), Arc(area_corners[2], area_corners[3]), Arc(area_corners[3], area_corners[0])] segment_boundaries = [Arc(segment_corners[0], segment_corners[1]), Arc(segment_corners[1], segment_corners[2]), Arc(segment_corners[2], segment_corners[3]), Arc(segment_corners[3], segment_corners[0])] angle1 = area_boundaries[0].angle(area_boundaries[1]) angle2 = segment_boundaries[0].angle(segment_boundaries[1]) if np.sign(angle1) != np.sign(angle2): segment_corners.reverse() segment_boundaries = [Arc(segment_corners[0], segment_corners[1]), Arc(segment_corners[1], segment_corners[2]), Arc(segment_corners[2], segment_corners[3]), Arc(segment_corners[3], segment_corners[0])] poly = [] boundaries = area_boundaries other_boundaries = segment_boundaries b__ = None for b__ in boundaries: if point_inside(b__.start, segment_corners): poly.append(b__.start) break else: inter = get_first_intersection(b__, other_boundaries) if inter is not None: poly.append(inter) break if len(poly) == 0: return None while len(poly) < 2 or poly[0] != poly[-1]: inter, b2_ = get_next_intersection(poly[-1], b__, other_boundaries) if inter is None: poly.append(b__.end) idx = (boundaries.index(b__) + 1) % len(boundaries) b__ = boundaries[idx] else: poly.append(inter) b__ = b2_ boundaries, other_boundaries = other_boundaries, boundaries return poly[:-1] R = 1 def get_area(corners): """Get the area of the convex area defined by *corners*. """ c1_ = corners[0] area = 0 for idx in range(1, len(corners) - 1): b1_ = Arc(c1_, corners[idx]) b2_ = Arc(c1_, corners[idx + 1]) b3_ = Arc(corners[idx], corners[idx + 1]) e__ = (abs(b1_.angle(b2_)) + abs(b2_.angle(b3_)) + abs(b3_.angle(b1_))) area += R ** 2 * e__ - math.pi return area def overlap_rate(swath_corners, area_corners): """Get how much a swath overlaps an area. """ area_area = get_area(area_corners) inter_area = get_area(polygon(area_corners, swath_corners)) return inter_area / area_area def min_distances(area_corners, segment_corners): """Min distances between each corner of *area_corners* and *segment_corners*. """ dists = np.ones(4) * np.infty for i, ic_ in enumerate(area_corners): for jc_ in segment_corners: dist = ic_.distance(jc_) if dists[i] > dist: dists[i] = dist return dists def should_wait(area_corners, segment_corners, previous_segment_corners): """Are the newest cornest still inside the area ? is the last segment boundary overlapping any boundary of the area ? In this case we should wait for the next segment to arrive. """ dists = min_distances(segment_corners, previous_segment_corners) indices = np.argsort(dists) new_corners = np.array(segment_corners)[indices[2:]] if len(new_corners) != 2: raise ValueError("More than 2 corners differ from previous segment...") new_arc = Arc(new_corners[0], new_corners[1]) for i in new_corners: if point_inside(i, area_corners): return True area_arc1 = Arc(area_corners[0], area_corners[1]) area_arc2 = Arc(area_corners[1], area_corners[2]) area_arc3 = Arc(area_corners[2], area_corners[3]) area_arc4 = Arc(area_corners[3], area_corners[0]) for i in (area_arc1, area_arc2, area_arc3, area_arc4): if i.intersects(new_arc): return True return False import unittest class TestSphereGeometry(unittest.TestCase): """Testing sphere geometry from this module. """ def test_angle(self): """Testing the angle value between two arcs. """ base = 0 p0_ = Coordinate(base, base) p1_ = Coordinate(base + 1, base) p2_ = Coordinate(base, base + 1) p3_ = Coordinate(base - 1, base) p4_ = Coordinate(base, base - 1) arc1 = Arc(p0_, p1_) arc2 = Arc(p0_, p2_) arc3 = Arc(p0_, p3_) arc4 = Arc(p0_, p4_) self.assertAlmostEqual(arc1.angle(arc2), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc2.angle(arc3), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc3.angle(arc4), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc4.angle(arc1), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc1.angle(arc4), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc4.angle(arc3), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc3.angle(arc2), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc2.angle(arc1), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc1.angle(arc3), math.pi, msg="this should be pi") self.assertAlmostEqual(arc3.angle(arc1), math.pi, msg="this should be pi") self.assertAlmostEqual(arc2.angle(arc4), math.pi, msg="this should be pi") self.assertAlmostEqual(arc4.angle(arc2), math.pi, msg="this should be pi") p5_ = Coordinate(base + 1, base + 1) p6_ = Coordinate(base - 1, base + 1) p7_ = Coordinate(base - 1, base - 1) p8_ = Coordinate(base + 1, base - 1) arc5 = Arc(p0_, p5_) arc6 = Arc(p0_, p6_) arc7 = Arc(p0_, p7_) arc8 = Arc(p0_, p8_) self.assertAlmostEqual(arc1.angle(arc5), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc5.angle(arc2), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc2.angle(arc6), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc6.angle(arc3), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc3.angle(arc7), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc7.angle(arc4), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc4.angle(arc8), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc8.angle(arc1), math.pi / 4, 3, msg="this should be pi/4") self.assertAlmostEqual(arc1.angle(arc6), 3 * math.pi / 4, 3, msg="this should be 3pi/4") c0_ = Coordinate(0, 180) c1_ = Coordinate(1, 180) c2_ = Coordinate(0, -179) c3_ = Coordinate(-1, -180) c4_ = Coordinate(0, 179) arc1 = Arc(c0_, c1_) arc2 = Arc(c0_, c2_) arc3 = Arc(c0_, c3_) arc4 = Arc(c0_, c4_) self.assertAlmostEqual(arc1.angle(arc2), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc2.angle(arc3), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc3.angle(arc4), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc4.angle(arc1), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc1.angle(arc4), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc4.angle(arc3), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc3.angle(arc2), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc2.angle(arc1), -math.pi / 2, msg="this should be -pi/2") # case of the north pole c0_ = Coordinate(90, 0) c1_ = Coordinate(89, 0) c2_ = Coordinate(89, -90) c3_ = Coordinate(89, 180) c4_ = Coordinate(89, 90) arc1 = Arc(c0_, c1_) arc2 = Arc(c0_, c2_) arc3 = Arc(c0_, c3_) arc4 = Arc(c0_, c4_) self.assertAlmostEqual(arc1.angle(arc2), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc2.angle(arc3), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc3.angle(arc4), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc4.angle(arc1), math.pi / 2, msg="this should be pi/2") self.assertAlmostEqual(arc1.angle(arc4), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc4.angle(arc3), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc3.angle(arc2), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(arc2.angle(arc1), -math.pi / 2, msg="this should be -pi/2") self.assertAlmostEqual(Arc(c1_, c2_).angle(arc1), math.pi/4, 3, msg="this should be pi/4") self.assertAlmostEqual(Arc(c4_, c3_).angle(arc4), -math.pi/4, 3, msg="this should be -pi/4") self.assertAlmostEqual(Arc(c1_, c4_).angle(arc1), -math.pi/4, 3, msg="this should be -pi/4") def test_inside(self): """Testing if a point is inside for other points. """ c1_ = Coordinate(-11, -11) c2_ = Coordinate(11, -11) c3_ = Coordinate(11, 11) c4_ = Coordinate(-11, 11) corners = [c1_, c2_, c3_, c4_] point = Coordinate(0, 0) self.assertTrue(point_inside(point, corners)) point = Coordinate(0, 12) self.assertFalse(point_inside(point, corners)) c1_ = Coordinate(-1, 180) c2_ = Coordinate(1, 179) c3_ = Coordinate(1, -179) c4_ = Coordinate(-1, -179) corners = [c1_, c2_, c3_, c4_] point = Coordinate(0, 180) self.assertTrue(point_inside(point, corners)) point = Coordinate(12, 180) self.assertFalse(point_inside(point, corners)) point = Coordinate(-12, 180) self.assertFalse(point_inside(point, corners)) point = Coordinate(0, 192) self.assertFalse(point_inside(point, corners)) point = Coordinate(0, -192) self.assertFalse(point_inside(point, corners)) # case of the north pole c1_ = Coordinate(89, 0) c2_ = Coordinate(89, 90) c3_ = Coordinate(89, 180) c4_ = Coordinate(89, -90) corners = [c1_, c2_, c3_, c4_] point = Coordinate(90, 90) self.assertTrue(point_inside(point, corners)) def test_intersects(self): """Test if two arcs intersect. """ p0_ = Coordinate(0, 0) p1_ = Coordinate(1, 0) p2_ = Coordinate(0, 1) p3_ = Coordinate(-1, 0) p4_ = Coordinate(0, -1) p5_ = Coordinate(1, 1) p6_ = Coordinate(-1, 1) arc13 = Arc(p1_, p3_) arc24 = Arc(p2_, p4_) arc32 = Arc(p3_, p2_) arc41 = Arc(p4_, p1_) arc40 = Arc(p4_, p0_) arc56 = Arc(p5_, p6_) arc45 = Arc(p4_, p5_) arc02 = Arc(p0_, p2_) arc35 = Arc(p3_, p5_) self.assertTrue(arc13.intersects(arc24)) self.assertFalse(arc32.intersects(arc41)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc45.intersects(arc02)) self.assertTrue(arc35.intersects(arc24)) p0_ = Coordinate(0, 180) p1_ = Coordinate(1, 180) p2_ = Coordinate(0, -179) p3_ = Coordinate(-1, -180) p4_ = Coordinate(0, 179) p5_ = Coordinate(1, -179) p6_ = Coordinate(-1, -179) arc13 = Arc(p1_, p3_) arc24 = Arc(p2_, p4_) arc32 = Arc(p3_, p2_) arc41 = Arc(p4_, p1_) arc40 = Arc(p4_, p0_) arc56 = Arc(p5_, p6_) arc45 = Arc(p4_, p5_) arc02 = Arc(p0_, p2_) arc35 = Arc(p3_, p5_) self.assertTrue(arc13.intersects(arc24)) self.assertFalse(arc32.intersects(arc41)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc45.intersects(arc02)) self.assertTrue(arc35.intersects(arc24)) # case of the north pole p0_ = Coordinate(90, 0) p1_ = Coordinate(89, 0) p2_ = Coordinate(89, 90) p3_ = Coordinate(89, 180) p4_ = Coordinate(89, -90) p5_ = Coordinate(89, 45) p6_ = Coordinate(89, 135) arc13 = Arc(p1_, p3_) arc24 = Arc(p2_, p4_) arc32 = Arc(p3_, p2_) arc41 = Arc(p4_, p1_) arc40 = Arc(p4_, p0_) arc56 = Arc(p5_, p6_) arc45 = Arc(p4_, p5_) arc02 = Arc(p0_, p2_) arc35 = Arc(p3_, p5_) self.assertTrue(arc13.intersects(arc24)) self.assertFalse(arc32.intersects(arc41)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc56.intersects(arc40)) self.assertFalse(arc45.intersects(arc02)) self.assertTrue(arc35.intersects(arc24)) def test_overlaps(self): """Test if two areas overlap. """ p1_ = Coordinate(89, 0) p2_ = Coordinate(89, 90) p3_ = Coordinate(89, 180) p4_ = Coordinate(89, -90) p5_ = Coordinate(89, 45) p6_ = Coordinate(89, 135) p7_ = Coordinate(89, -135) p8_ = Coordinate(89, -45) self.assertTrue(overlaps([p1_, p2_, p3_, p4_], [p5_, p6_, p7_, p8_])) self.assertFalse(overlaps([p1_, p5_, p2_, p6_], [p3_, p7_, p4_, p8_])) p1_ = Coordinate(1, -1) p2_ = Coordinate(1, 1) p3_ = Coordinate(-1, 1) p4_ = Coordinate(-1, -1) p5_ = Coordinate(0, 0) p6_ = Coordinate(0, 2) p7_ = Coordinate(2, 2) p8_ = Coordinate(2, 0) self.assertTrue(overlaps([p1_, p2_, p3_, p4_], [p5_, p6_, p7_, p8_])) self.assertFalse(overlaps([p1_, p8_, p5_, p4_], [p2_, p3_, p6_, p7_])) def test_overlap_rate(self): """Test how much two areas overlap. """ p1_ = Coordinate(1, -1) p2_ = Coordinate(1, 1) p3_ = Coordinate(-1, 1) p4_ = Coordinate(-1, -1) p5_ = Coordinate(0, 0) p6_ = Coordinate(0, 2) p7_ = Coordinate(2, 2) p8_ = Coordinate(2, 0) self.assertAlmostEqual(overlap_rate([p1_, p2_, p3_, p4_], [p5_, p6_, p7_, p8_]), 0.25, 3) c1_ = [(60.5944, 82.829699999999974), (52.859999999999999, 36.888300000000001), (66.7547, 2.8773), (80.395899999999997, 98.145499999999984)] c2_ = [(62.953206630716465, 7.8098183315148422), (62.953206630716465, 26.189349044600252), (53.301561187195546, 26.189349044600252), (53.301561187195546, 7.8098183315148422)] cor1 = [Coordinate(t[0], t[1]) for t in c1_] cor2 = [Coordinate(t[0], t[1]) for t in c2_] self.assertAlmostEqual(overlap_rate(cor1, cor2), 0.07, 2) c1_ = [(60.5944, 82.829699999999974), (52.859999999999999, 36.888300000000001), (66.7547, 2.8773), (80.395899999999997, 98.145499999999984)] c2_ = [(65.98228561983025, 12.108984194981202), (65.98228561983025, 30.490647126520301), (57.304862819933433, 30.490647126520301), (57.304862819933433, 12.108984194981202)] cor1 = [Coordinate(t[0], t[1]) for t in c1_] cor2 = [Coordinate(t[0], t[1]) for t in c2_] self.assertAlmostEqual(overlap_rate(cor1, cor2), 0.5, 2) if __name__ == '__main__': unittest.main()