First bunch of files
[ais.git] / bin / area.py
1 #!/usr/bin/env python
2 # -*- encoding: utf-8
3
4 __all__ = [ 'Area', 'load_area_from_kml_polygon' ]
5
6 import sys
7
8 class Area:
9     """
10     That class defines an area (on the Earth)
11     It provides testing whether a point is inside or not
12     """
13     def __init__(self, points=[]):
14         self.points = []
15         for p in points:
16             self.addpoint(p)
17     def addpoint(self, point):
18         self.points.append(point)
19         if len(self.points)==1:
20             self.min = point
21             self.max = point
22             return
23         if point[0] < self.min[0]:
24             self.min = (point[0], self.min[1])
25         elif point[0] > self.max[0]:
26             self.max = (point[0], self.max[1])
27         if point[1] < self.min[1]:
28             self.min = (self.min[0], point[1])
29         elif point[1] > self.max[1]:
30             self.max = (self.max[0], point[1])
31
32     def check(self):
33         for point in self.points:
34             if not self.contains(point):
35                 return False
36         return True
37
38     def contains(self, point):
39         if not self.points:
40             return False
41         # first test the bounding box
42         #if point[0] < self.min[0] \
43         # or point[0] > self.max[0] \
44         # or point[1] < self.min[1] \
45         # or point[1] > self.max[1] :
46         #    return False
47         for i in range(len(self.points)):
48             p1 = self.points[i]
49             x1, y1 = p1
50             p2 = self.points[(i+1)%len(self.points)]
51             x2, y2 = p2
52             xa = point[0] - x1
53             ya = point[1] - y1
54             xb = x2 - x1
55             yb = y2 - y1
56             if xa * yb < xb * ya:
57                 return False
58         return True
59
60 def load_area_from_kml_polygon(filename):
61     file = open(filename)
62     coordinates_lines = [ line for line in file.readlines() if '</coordinates>' in line ]
63     if len(coordinates_lines) != 1:
64         print >> sys.stderr, 'There should be exactly one line with coordinates in', filename
65         sys.exit(1)
66     coordinates = coordinates_lines[0].replace('</coordinates>', '').replace('\n', '').replace('\r', '')
67     coordinates = [ xyz for xyz in coordinates.split(' ') if xyz ]
68     if coordinates[0] != coordinates[-1]:
69         print >> sys.stderr, 'First and last coordinates of', filename, 'should be the same'
70         print >> sys.stderr, coordinates[0]
71         print >> sys.stderr, coordinates[-1]
72         sys.exit(1)
73     assert len(coordinates)>3, 'polygon should have 3 edges minimum'
74     
75     area = Area()
76     for xyz in coordinates[0:-1]:
77         x,y,z = xyz.split(',')
78         area.addpoint((float(y),float(x)))
79     assert area.check(), 'Polygon should be counter-clockwise and convex.'
80     return area
81
82
83 #if __name__ == '__main__':
84 # counter clock-wise : Positive
85 #pelagos = Area([
86 #    (42.91, 12.5),
87 #    (45.3612930132714, 10.01843703552244),
88 #    (43.6,5.5),
89 #    (40.57,8.6)
90 #    ])
91 #for p in [
92 #    (42,9),
93 #    (41,5),
94 #    (40,12),
95 #    (45,13),
96 #    (45,7),
97 #    ]:
98 #    print "testing", p
99 #    if pelagos.contains(p):
100 #        print "inside"
101 #    else:
102 #        print"outside"
103
104